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A STABILITY ANALYSIS OF CYLINDRICAL PANELS 
USING A FINITE ELEMENT FORMULATION 

by 

Richard E. Snyder 
ABSTRACT 

A cylindrical finite element suitable for the linear stability 
analysis of cylindrical shells is developed. Energy principles and vari- 
ational methods lead to a problem formulation which lends itself to 
physical interpretations of the governing matrices of the finite element. 
By properly grouping the terms which result from taking the second varia- 
tion of the potential energy of the element > it is possible to identify 
three distinct types of matrices. The first matrix is the conventional 
stiffness matrix; the second is an "initial stress" stiffness matrix; 
and the third is an "initial displacement" stiffness matrix. With the 
assumption of linearity, the buckling problem is stated in terms of the 
classical linear real eigenvalue equation. This problem formulation was 
programmed on the CDC 6600 series computer. The computer program is used 
to analyze the buckling of a variety of structures. Columns, arches, flat 
plates and curved panels with and without cutouts are considered. Com- 
parisons are made between closed form solutions and the results of the 
present analysis to establish confidence in the techniques used. Curved 
panels with cutouts of varying size are analyzed for buckling. The in- 
fluence of curvature and cutout size on the prebuckling deformations in a 
curved panel are studied and found to be significant. The prebuckling 
deformations are shown to have a significant influence on the buckling 
strength of curved panels with cutouts. 
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CHAPTER X. 


INTRODUCTION 

Plate and shell stability analysis and its application to practical 
engineering structures has been the subject of extensive structural re- 
search efforts. Despite this, many problems still remain in the accurate 
prediction of the buckling mode of failure of many types of practical 
plate and shell configurations, particularly those with holes or cutouts. 
For example, reliable procedures for the analysis of a flat plate with a 
large cutout are limited. The finite element method provides a means of 
analyzing problems such as these. The development of this method has 
been prompted by its ability to model complex geometries and loading 
conditions. To date, major developments in finite element methods have 
concentrated on stress and vibration analysis with only limited attention 
given to stability. The objective of this dissertation is to report the 
results of an extension of the finite element method to problems of thin 
shell instability. The extension will concentrate on the development 
and use of elements which, thus far, have seen only limited application 
in finite element technology. Primary emphasis will be placed on the 
buckling analysis of curved panels with rectangular cutouts. However, 
the method will also be applied to columns, flat plates and arches. 

The field of shell stability analysis is one of the most extensively 
investigated areas of classical mechanics. The list of references 
associated with classical shell stability is very formidable. References 
to classical shell buckling investigations are cited herein only to the 
extent that they relate to special aspects of the finite element approach. 
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Formalized finite element methods generally began with the presenta- 
tion of the "direct stiffness" method by Turner, Clough, Martin, and 
Topp (ref. 1). Since then, finite elements have been used in the stress 
analysis of a wide variety of structures, including frames, arches, 
shells, and solids. In recent years, attention has been turned to the 
application of finite element technology to the stability of structures. 
Expansion of the stiffness method to handle nonlinear, large deflection 
problems was first presented in 1960 by Turner, Dill, Martin and Melosh 
(ref. 2). In 1962, Turner, et.al., (ref. 3) enlarged on the nonlinear 
finite element technique by presenting an eigenvalue procedure to deter- 
mine the buckling of columns. Gallagher and Padlog (ref. 4) independent- 
ly derived a stability coefficient matrix to predict column buckling. 
Summaries of the current state of the art of stability predictions using 
beam-column elements are contained in references 5, 6 and 7. In these 
references, the stability problem is framed in terms of the conventional 
stiffness matrix and what is termed the "geometrical" or "initial stress" 
stiffness matrix. The term "initial stress" stiffness matrix will be 
adopted in this work. The terminology applied to this new matrix re- 
flects its dependence on the initial state of stress and undeformed 
geometry of the element. 

In reference 8, Gallagher, et.al., extended finite element stability 
methods to flat triangular elements. The explicit formulation of the 
initial stress matrix for a rectangular plate in bending is presented . 
by Kapur and Hartz in reference 9. The stability of doubly curved shells 
of revolution subjected to axisymmetric loading, using the finite element 
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method, was investigated by Navaratna, Pian and Witmer (ref. 10). The 
stability of cylindrical shells using curved finite elements was studied 
by Bogner, Fox and Schmit (ref. 11). In reference 11, the problem is 
formulated from the standpoint of direct minimization of the total 
potential energy as opposed to the development of identifiable stiffness 
and initial stress stiffness matrices. In addition, a large number of 
degrees of freedom (i.e.', 48) are used for each element. 

There is limited literature published on the buckling of curved 
cylindrical panels. Classical analyses of curved panels typically con- 
sider the case of an infinate aspect ratio (ref. 12 and 13). Gerard and 
Becker (ref. 14) directly consider "very wide" and "very narrow" curved 
panels and then fair a curve between those results to cover panels of 
intermediate dimensions. The importance of boundarv conditions in the 
determination of the buckling of curved panels is established by 
Rehfield and Hallaur (ref. 15). 

The effect of cutouts in cylinders has been the subject of several 
experimental investigations (ref. 16* 17 and 18). Brogan and Almroth 
(ref. 17) applied a finite difference approximation to the governing 
equations and obtained reasonable agreement with experiments. No re- 
ferences dealing with the bifurcation buckling of curved panels with 
cutouts were found, and indicates a need for data on such problems. 
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CHAPTER II. 

APPLICATION OF ENERGY PRINCIPLES AND 
VARIATIONAL METHODS TO LINEAR STABILITY ANALYSIS 

II. -1 Basic Principles 

The principle of minimum potential energy establishes that for 
equilibrium the total potential energy, tt , for a system must be extremal 
or stationary (ref. 19 and 1 20) . Thus for equilibrium, the first varia- 
tion of the potential energy vanishes. 

(Sit = 0 (II-l) 

where 6 is the variational symbol. 

The stability of the equilibrium state can be investigated by ex- 
amining the second variation of the potential energy. An equilibrium 
state is stable if every neighboring state has a larger potential energy 

In other words, an equilibrium state is stable if, in addition to satis- 

o 

fying equation II-l it also satisfies the condition 5 i '7T>0. Conversely, 

_ 2 

equilibrium is unstable if 6 tt < 0. Therefore, the infinitesimal 
stability limit as used herein corresponds to the case of 

S 2 tt = 0 (II-2) 

This principle is well-known for continium problems and has been 
applied to approximations based on finite element methods in references 
4, 5, 21, 22, 23 and 24. A formulation of the second variation of the 
potential energy that lends itself conveniently to the numerical approxi- 
mations found in the finite element method will be presented. In this 
chapter, these concepts will be applied to a beam-column to illustrate 
the approach. In Chapter III, the same methods will be applied to the 
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more complex curved panel problem. 


II. -2 The Beam Column 

The beam-column finite element shown In Figure II-l extends be 


tween node points 1 and 2, and has a cross sectional area. A; a 
length £; a moment of inertia, I; and a modulus of elasticity, E. 
The forces acting at each end of the beam— column are shown in their 
positive directions in Figure II-l. The displacement in the Z direc- 
tion, w, and the displacement in the X direction, u, are also 
shown in their positive directions. The total potential energy of the 
finite element system is 


H & 2 

7T = y / EAe'aX + / w, xx dX- -NjU " M x 0 l ~ N 2 U 2 " V 2 W 2 ~ M 2 6 2 

z 0 2 


0 


where (reference 19) 


e = c + — 
2 


(II- 3) 
(II-4) 


and 


e = u, x (II-5) 

<t> = w, (II-6) 

r *x 

Here, e is the nonlinear middle surface strain composed of the linear 
strain £ and the rotation <j) . The first integral on the right hand 
side of equation II-3 is the membrane strain energy and the second in- 
tegral is the bending strain energy. The remaining terms represent the 
potential energy of the external forces. 



- 6 - 


Z, w 



E = Modulus of elasticity 
A = Area 

I = Area moment of inertia 


Figure II- 1. The Beam Column 
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Substituting equation II-4 into II-3 and rearranging leads to 

L L L t L 2 

IT = — / EAe 2 dX + I / EA<jr(e+ i^ 2 )dX - I / EA* dX + i / EIw, dX 

20 2 0 2 204 2 0 

— N^u^ — V^w^- ^2^*2 ^2^2 — ^2^2 (1^ — 2) 


The axial force and bending moment may be written as 

N = AEe - AE (e + -| 0 2 ) (Il“8) 

M = El w } (II-9 ) 

AA 

Substituting II-8 and II-9 into II-7, 

£ £ £ , £ 
tt = I / EAe aX + 1 f N^ 2 dX - ~ f EAcJ> dX + “ / Mw. dX 
2 0 2 0 2 0 4" 2 0 

— N-^u^ — V^w^ — M-^0^ ~ N 2 u 2 — V 2 W 2 — M 2 0 2 (11-10) 


After giving the system a virtual displacement, the total potential 
energy becomes 


tt + Att = ~ f EA(e+6e) 2 dX + — / 
2 0 2 0 
n 

(N+AN) (<£>+6<j>) 2 d X - - f M(fh5(^) 4 dX 

2 0 4 

1 * 

+ j / (M+AM) (w, xx +6w, xx )dX - N x u x - V jv^ - M 1 0 1 - N 2 u 2 - V 2 w 2 

- M 2 0 2 - N x 5u x - V x 6 Wl - 

M l 60 l“ N 2 6 u 2 - V 2 6w 2 - M 2 6G 2 (11-11) 

where 

r 2 1 


an =ea jAe + 

(11-12) 

AM = El <5 w, 

XX 

(11-13) 


Substituting 11-12 and 11-13 into 11-11, 
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TT + ATT = i /eA(c 2 +2£«E + 6 2 e)dX + i / N(<f> 2 + i 06<f + 5 2 4>) dX 
2 0 2 0 2 

+ i (4> 2 + + fi 2 <t>)dX 

1 0 \ 2 / 

£ 

- i / EA (<j> 4 +4<t> 3 d<j> + 6tf> 2 6 2 4> + 4<|>6 3 <t> + 6 4 <(i)dX 
8 0 

i 2 1 2 i 2 

+ \ f Q Mw -xx dX + '| f Q M<Sw ’xx dX + f f Q EIw ’xx 5w . x x dX 


(11-14) 


+ — f El6 2 w, xx dX - N^u-^ - Vjw-j^ - - N 2 U 2 “ ^2 X/J 2 “ 

2 0 


- N-^du^- V-^6 w^ - M-^60^ - N 2 ^ U 2 “ ^ 2^ w 2 ” ^ 2^^ 2 
rr + At t may be expressed in the following form (ref. 19) 


1 9 

7T ~h ATT = TT + 67T + — 6 TT H- — — 

2! 


(11-15) 


By arranging equation 11-14 into the form of equation 11-15, 
the matrix form of the second variation is 


6 2 tt 


f [6c 64> <5w, xx ] 
0 



"100" 


000“ 


*°4> o" 


6c 

EA 

000 

+ 

0N0 

+ EA 

U 0 


64> 

- 

OOJ 
L AJ 

. 

000 


00 0 

L J 


6w >xx 


dx (11-16) 


Equation 11-16 may be stated more conveniently as 

<5 2 tt - / {Se} T [K°]{6c}dx + f {6e} T [K 1 ] {6e}dx + / (6e} T [K 2 H6e}dx 
0 0 0 


(11-17) 


where 


{5e}- 


6e 

6<$i 

6w, xx 


9 



Prior to introducing the numerical approximations involved in the 
finite element method, several general statements mat he made concerning 
the use of equation 11-17 for finite element type buckling analyses. 

0 Terms relating to the work of the external forces do not appear 
in the second variation of the potential energy. 

0 [K°] leads to the conventional stiffness matrix for a beam- 

column. (See reference 24 and Appendix A). 

0 [K-*-] is a function of the load in the element and leads to a 

matrix denoted herein as the "initial stress" stiffness matrix. 

° [K 2 ] is a function of rotational deformations and leads to a 

matrix denoted herein as the "initial displacement" stiffness matrix. 

° The finite element numerical approximations (i.e., displace- 
ment function) need not necessarily be the same when dealing with 
[K^], [K^] and [K 2 ]. However, if convergence is to be obtained in the 
limit, all of the numerical approximations must represent the essential 
character of the problem. 

0 A static solution to the problem must first be accomplished to 
establish values of N and which are required for the evaluation of 
[K 1 ] and [K 2 ]. 

The details of the steps required to develop the conventional, 
the initial stress and the initial displacement stiffness matrices for 
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both the beam-column and the arch are given in Appendix A. The matrix 
formulation of the strain-displacement relations and the assumed dis- 
placement function as well as intermediate matrix products are presented. 
The symbols used in Appendix A correspond to those used in this chapter. 

For the beam-column, the strains are related to the displacements, 
u and w, through a matrix of differential operators. 

(B = [D] {g} <n- 18 > 

A critical feature of any finite element development is the dis- 
placement function chosen to represent the deformation characteristics 
of the element in terms of the nodal displacements. A linear variation 
of u and a cubic variation of w, as is used in Appendix A, is fre- 
quently assumed (ref. 5). Such a set of assumptions stated in matrix 
form is 

{g} = [B] [r] {a} (11-19) 

The matrix [B] is a function of X and the [F] matrix is a function of 
the element geometry. {A} is the vector of nodal displacements. 
Substituting equation 11-19 into equation 11-18, 

{7} = (G) {a} (11-20) 

where 

[G] = [D] [B] [ r ] 
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where L 

[K° ] = / [G] T [K°] [G] dX 
e 0 

L 

[K X e ] = / [G1 T [Kl] [G] dX 
0 

L 

[K 2 ] = / [G] T [K 2 ] [G] dX 
e 0 

In general, A, E and I, are functions of the longitudinal 
coordinate, X. This is also true of the load, N, and the rotation, 

<j> . Before the integration of the terms in equation 11-21 may be 
carried out, an assumption must be made regarding the variation of these 
parameters along the length of the beam-column. The assumption that 
A, E and I are constants over the length of the element is commonly 

i 

made. This assumption will be adopted here. Further, N and <j> will 
also be assumed to be constants over the length of the element. The 
matrices [K® ], [K^ ] and [K 2 ] given in Appendix A are based 

on these assumptions. 

Applying the neutral stability condition, <5 z tt = 0, the buckling 
criterion for the, beam-column element is 

det [K° e ] + lK 1 e ] + [K 2 e ] = 0 (11-22) 

The eigenvalue involved in the solution of equation 11-22 is 
the ratio of the bifurcation load, N cr , to the applied load, N. 

The load and rotation at buckling are equal to the initial load and 
rotation times the eigenvalue. Because of the <t>^ term in [K , 
equation 11-22 has the form of a quadratic eigenvalue problem. 
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The nonlinear strain displacement relations used in this develop- 
ment assume the square of the rotations to be of the same order as the 

2 

strains, and the strain to be small compared to unity or e = 0 (<f> ) « 1 » 

Thus, a reasonable first approximation to the solution of equation 11-22 

for many problems is to assume the $ terms to be negligible compared 

to the 4) terms. With this approximation, equation 11-22 becomes a 

1 

linear eigenvalue problem and is 




r* — 


det 

[K° e ] + X 

[Kl e ] + [K 2 e ] 
— — 



where 


(11-23) 


X = N cr /N 

Ner = buckling load 

N = applied load 


The method of assembling element stiffness, initial stress stiff- 
ness and initial displacement stiffness matrices to represent a complete 
structure will be discussed in Chapter III. The method used in the 
determination of the eigenvalues is described in Chapter IV. 
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CHAPTER III. 

TH E CYLINDRICAL, THIN SHELL, FINITE ELEMENT 

The following development is an extension of a stiffness formula- 
tion by Gallagher (ref. 25) for the doubly curved shell element, shown 
in Figure III-l, to include elastic instability effects. The Gallagher 
shell element was selected as a basis for this work because it has been 
shown to give reliable results and because it is well documented. An 
expression for the strain energy of the doubly curved shell element will 
first be derived. Then, prior to introducing an assumed displacement 
function, the problem will be specialized to a singly curved cylindrical 
element. 

The cylindrical element has wide application in aerospace type 
structures. It can be used in the analysis of structures such as 
airplane fuselages, rocket motor cases, tanks, and interstage adaptors. 
Elliptical cross sections may be represented with cylindrical elements 
by allowing the radius of curvature to vary from element to element. 

III.-l Basic Assumptions 

In this development the shell material is assumed to be isotropic 
and to obey Hooke's Law. The neutral surface of the shell lies midway 
through the thickness. Applying the Kirchhof f-Love hypothesis, it is 
further assumed that: 

1. The displacements u, and v corresponding to the directions 
and ^2 (Figure III-l) respectively, are linear in the 
thickness direction, £ 3 . 
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2. All components of stress normal to the middle surface are 
negligible. 

3. The displacement normal to the middle surface is a function 
only of the middle surface coordinates. 

HI. -2 Formulation of the Potential Energy 

As was illustrated in Chapter II, the potential energy of the 
external forces applied to an element does not appear in the expression 
for the second variation of the potential energy. Hence, consideration 
will be directed only to the evaluation of the strain energy of the 
element . 

The geometry of the doubly curved, thin shell, finite element is 
depicted in Figure Ill-la. The middle surface of the element is defined 
by the curvilinear coordinate and £ 2 * T ^ e coordinate £3 is 

normal to the middle surface and completes the orthoginal right-handed 
system. The radii of curvature R x and R 2 , corresponding to the 
coordinate lines and £ 2 respectively, are constants. The equa- 

tion for the differential distance, ds, between two points on the 
middle surface is 

ds 2 = a ] 2 d ^ 2 + a 2 2 dS 2 2 • (III-D 

where a-^ and a 2 are the Lame parameters. The linear displacements 
u , v, and w, corresponding to the coordinate directions £2 

and £3 respectively, are as shown in Figure III- lb. 
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The nonlinear strain displacement relations for small strains 
were derived by Sanders in reference 26 and will be used in the formula- 
tion of the strain energy expression. The nonlinear middle surface ex- 
tens ional and in plane shear strains are; 


3i 2 

e l " e l + ~T” 

32 2 

e 2 “ £ 2 + ~J~ 


(III-2) 


e 12 88 >12 + frl &2 

where 


e l 


i + w 

oil R-i 


= i_ 8 V , 

2 a 2 8^ 2^ r 2 


y _ ^ . _1 5U 

12 “ a l <*€ 1 <*2 ?2 


D 1 8w 

6 1- + 


u 
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The expressions for the bending distortion are: 
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The expression for the strain energy for the doubly curved shell 
seleraent of Figure Ill-la, expressed as integrals over the middle surface 
area, is (ref. 19) 


U = 


Eh 


2(l-\)2) area 


[ e i 


2 o 

// | e, +e 2 + 2ve-,e„ + 


(1-v) 2 

12 T 2 12 
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(l-V) y P 
2 12 


a^d^d^ 

a l a 2 d ^l d ^2 


(III-4) 


The total strain energy is seen to be the sum of the membrane energy, 
given by the first integral, and the bending energy, given by the second 
integral. 

The stress components in the shell element are shown in Figure 
III-2a. The corresponding stress resultants and bending moments are 
shown in their positive directions in figures III-2b and III-2c. 

The following set of equations relate the stress resultants to membrane 
strains and the. bending moments to bending distortions (ref. 19). 


(e-^ + ve 2 ) 



(III-5) 
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By combining equations III-2 and III-5 with equation III 4, the follow- 


ing form of the strain energy expression is obtained: 


U = I // 
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HI. -3 The Second Variation of the Potential Energy 

The strain energy of the shell element after a virtual displace- 
ment may be expressed as (ref. 19) 


U + AU = U + 6U+ — 6^11 + High Order Terms (III-7) 

2! 

9 

in which 6U and 6 U are the first and second variations, respective- 
ly, of the strain energy. . Hence, following the technique used in 
Chapter II, the second variation of the strain energy will be determined 
by giving the system a virtual displacement and grouping the terms in 
the resulting strain energy expression in the form of equation III-7. 
Applying a virtual displacement to the terms in equation III-7 produces 
the following equation: 


U + AU = i J/ a | (Nj+ANi) 
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Evaluating the incremental stress resultants and moments: 
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Comparing this with the first of equations III-5, it is determined that: 
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Similarily 
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After substituting equations III —9 through III-14 into equation III-8 
and after considerable manipulation, the second variation of the 
potential energy, which is equal to the second variation of the strain 
energy, may be identified and written in matrix form as: 


21 


5 2 tt „ ff {67} t — [^HSUa^d^dSa + //{«e} T [K 1 H«e}<x 1 oi 2 dC 1 d5 2 
1-V 2 area 

area 

+ //{ 6 e} — [K 2 ] 012^1^2 

1-V 2 

area 


where 


{67} t = ] Jc 1 63i 6X X 6^12 6£ 2 6 &2 ^ x 2 6x 12j 


[K°] - 


m - 


[K 2 ] 


6Ej 

1 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 


63 


0 

0 

0 

0 

0 

Nl 

.0 

0 

0 

N 

0 

0 


6X-l &y±2 ^ 


12 


0 3x 

3i ^i 2+ (^) B 2 i 

0 0 

» (i=y2 

0 v3i 

v3 2 &1&2 

I 0 o 

Lo o 


0 

0 

h 2 ■ 
12 
0 

0 

0 

Vh 2 

12 

0 

0 

0 

0 

0 

0 

JO- 

0 

0 


0 

0 

0 

1-V 

2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


63^ 

t 

0 

0 

0 


0 

N 12 

0 

0 

0 

N2 

0 

0 

v32 


° ("MV 2 V ^- 


6X2 


0 

0 

Vh 2 

12 

0 


0 (¥H 


(¥) Sl 

$2 

S2 2 +(¥) Bi: 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


(III-15) 


fix 


12 

0 

0 

0 


0 0 

0 0 

h 2 0 

12 h 2 (l-v) 
0 


24 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


1 



22 


The first term in Equation III-15 leads to the conventional stiffness 
matrix; the second term leads to the initial stress stiffness matrix; 
and the third term leads to the initial displacement stiffness matrix. 


III. -4 .The Stiffness Matrix 

The formulation of a stiffness matrix for the doubly curved shell 
element shown in Figure Ill-la is described in detail in reference 25 
and will not be repeated here. The stiffness matrix reported in refer- 
ence 25 specialized to the case of the element shown in Figure III-3a, 
will be used. The relationship of the coordinate axes of the cylindri- 
cal element to those of the doubly curved element is: 


= X (in- 16 ) 

C 2 - e 

5 3 - z 


The radii of curvature for the cylindrical element are: 


(III-17) 


R 2 - R 


The Lame parameters for the cylindrical element are: 

a i ■ i 

(III-18) 

c*2 = R 

Figures III-3a and III-3b depict the positive directions for the linear 
and angular displacement in the cylindrical element. 

Because of the importance of the assumed displacement function in 
the development of any finite element stiffness matrix; it is worth- 
while to briefly discuss the displacement functions used by Gallagher 


a. Element Geometry 


Z 



h. Linear Displacements c, Augu l&r Displacements 


Figure 1TI-3. Geomertry of the cylindrical element and associated 

displacements a 
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for the doubly curved shell element. The displacement functions used 
by Gallagher in reference 25 are: 

u = JL [ (X-a) (R0-b)u, - X(R0-b)u 2 + XR0u 3 - R0(X-a)u4] 
ab 

v = -i. [(X-a)(R0-b)v 1 - X(R0-b)v 2 + XR0v 3 - R0(X-a)v 4 ] 
ab x 

... ■ 1 I (a 3 + 2X 3 -3aX 2 )[b 3 +2(R6)3-3b(Re) 2 ]M 1 +(3aX 2 -2X 3 )[b3+2(Re) 3 

a 3b3 \ 

-3bR0) 2 ] w 2 

+ (3aX 2 -2X 3 ) [3b(R8) 2 -2(Re) 3 ]«3+(a 3 +2X 3 -3aX 2 ) [3b(R8) 2 -2(R9) 3 ]w 4 

+ aX(X-a) 2 [b 3 -2(Re) 3 -3b(R9) 2 J<t'x 1 +art 3 -aX 3 )[b 3 +2(Re) 3 -3b(Re) 2 ]<(i X2 

+ a(X 3 -aX 2 ) [3b(R0) 2 -2(R8) 3 ]<t x ^+a(X-a) 2 X[3b(R6) 2 -2(R8) 3 ]ij>^ 

+ b(a 3 +2X 3 -3aX 2 )R0[(R8)-b] 2 <{i0 -b(3aX 2 -2X 3 )R9[(R0)-b] 2 (}>Q 2 

+ b(3aX 2 -2X 3 ) [(R9)3-b(R6) 2 ]c)>e 3 +b(a 3 +2X 3 -3aX 2 )((R0) 3 -b(R9) 2 ]iti e4 

+ abXR0(X-a) 2 [ (R6)-b ] 2 <l>x0 +abXR0(X 2 -aX) [ (R0) -b ] 2 <t> x0 

l ^ 

+ abXR0(X 2 -aX)[(R0) 2 -bR0]<(i x9 +abXR0(X-a) 2 t(Re) 2 -bR6]i)) x04 j 

where, as shown in Figure III-3,a and b are the element lengths in the 

meridinal and hoop directions respectively; u^, and w^ (i=l,2,3,4) 

are the linear displacements at the i th corner of the element and; 

a A and 6 are the angular displacements at the i Ln 

t X £ 0 ± XOi 


corner of the element. 
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Zienkiewicz and Cheung in reference 27 list the desirable condi- 
tions to be met by a displacement function chosen to represent element 
behavior. These conditions are that the displacement function must 
properly account for rigid body motion and constant strain rates, 
and must satisfy inter-element boundaries, ’The above displacement 
functions meet these conditions in the case of a flat plate but fail to 
do so in the case of the curved element. 

Previous studies (ref. 28 and 29) indicate that the violation of 
the above conditions does not prevent convergence to the classical 
solution and does not significantly reduce accuracy for refined idealiza- 
tions. Indeed, Gallagher demonstrates the adequacy of his formulation 
by showing excellent corelations with known closed form solutions to 
several shell problems. 

Table III-l shows the organization of the terms in the conven- 
tional, initial stress and initial displacement stiffness matrices. 

The explicit statement of the terms in the element stiffness matrix, 

[K° 3 , obtained by Gallagher and specialized to the cylindrical element 
is given in Table III-2. 
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Table III-2 (Continued) 


k i 7 :x 

= k 

23,1 

= k 13 

,11 

0 

0 

0 

k 18,l 

= k l9 ,6 

" k 24 

,7 


0 




0 

= Dm 

b 

k 19 ,1 

= k 13,7 

6a 


“20 


= k H 


b VD 


m 


= ik 


0 

5,1 


6b 


0 0 0 , ° _ t 0 
r ^ V n — k- J> 1 “ kn A 1 / k( 


,1 " k 7,2 " k 13,8 k l9,l4 k 8,l 


k 22 ,1 

= k 19 , 4 

” k 16 ,7 

k 13,10 

“ k 16,l 




0 

0 

0 

0 

0 




k 24 ,1 

= k 

13,6 

= k 18 , 7 

“ k 19 , 12 

_k i8,; 









_ 

1 


0 

0 

0 

0 

f a 


Db 

1-v b 

k 2 ,2 

= k 8,8 : 

= k 14, 14 

= k 20 , 20 = 


i 

_ 

+ — - 

D m R 2 

6 a 

0 

0 

-7aD m 

(2-3v) 





k 3,2 

" k 9,g = 

+ 

40R 

- Dk 
2aR D 






1+4 

L 


0 

: 4 ,2 

0 

'5,2 

0 

"6,2 

0 

c 8 ,2 

0 

*9,2 


0 

0 

0 

a 

= k 22,2 

“ k 14 ,10 = 

= k 16, 14 

40R 

0 

0 

0 

7abD m 

= k ll,8 

= k 17 , 14 : 

= k 23, 20 

240R 




'l 

0 

0 

0 

a b ° m 
m 

= k 18,8 

= k 14 ,12 : 

= k 24 ,20 

240R 

0 

f a 

Db “| 

b 

= k„„ . 

= D m ( — 

1+ ~rT , 

- (1-v) - 

20,4 

m |6b 

L v . 

a 

0 

3 a Dm 

(2-3v)D b 

= k 8,3 

40R 

2aR 



20bR 


1+4 


% 


Db 

DmR 2 jJ 


28 


Table III - 2 (Continued) 
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I I 1-5 The Initial Stress a nd Initial Displacement Stiffness Matrices 

The steps involved in the development of the initial stress and 
initial displacement stiffness matrices for a cylindrical element will 
be outlined in this section. The detailed statement of the significant 
intermediate matrices involved in this development is in Appendix B. 

The matrix symbols used in this section agree with those in Appendix B. 

The desired strain-displacement relations for a cylindrical shell 
are obtained from equations III-2 and III-3 by the direct substitution 
of equations III-16, III— 17 and III— 18. Written in matrix notation, the 
equations relating strains to linear displacements are: 

{el = [A ] [D ] ( g } (III-2Q) 

where 

T 

{g } = U V wj 

The terms in the [A] matrix are all constants and [D] is a matrix of 
differential operators. 

One of the major assumptions involved in the development of the 
initial stress and initial displacement stiffness matrices is the form 
of the displacement functions to be used. Three displacement components, 
u, v and w must be characterised. The characterization of the membrane 
displacements, u and v, is based on the simple assumption of linear edge 
displacements. For a flat plate, this assumption insures compatibility 
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of displacements along lines bounding the elements. Making use of the 
approach taken in references 30 and 31, the displacement functions for 
u and v are: 


u 


v 


~ [(X-a) (R0-b ) Ul -X (R8-b) u 2 + XR0 u ? -R 6 (X-a) u 4 J 
^[(X-a) (Re-b)v 1 -X(R0-b)v 2 + XRdv^Re (X-a) v 4 J 


( 111 - 21 ) 

(III-22) 


In reference 6, Martin discusses the relative merits of linear 
versus cubic displacement functions for the case of a beam-column. He 
demonstrates that a linear function representing the normal displacement, 
w, of the beam-column is. the simplest, nontrivial polynomical form con- 
sistent with the problem. The beam-column stability problem is 
formulated using a cubic displacement function for w to derive the 
conventional stiffness matrix and a linear displacement function for w 
to derive the initial stress stiffness matrix. This is effectively a 
superposition of a tension-compression member and a beam; with no inter- 
action between the two. in the case of the beam, this has led to satis- 
factory results. In an analogous manner, the displacement function 
chosen for w in the case of the cylindrical element is: 


w - [(X-a) (Re-b)w 1 -X(R0-b)w 2 +XROw 3 -R0(X-a)w 4 J (III-23) 


The displacement functions, in matrix notation are: 


{g} = [B]{A} 


(III-24) 
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Where : 

{A} T -L“l “2 u 2 u 4 V 1 v 2 v 3 v 4 W 1 w 2 w 3 w 4_j 

and [B] is given in Appendix B. 

Substituting equation III-24 into equation III-20 gives: 

{e} = [G ] { A} (III-25) 

with [ G ] being given in Appendix B. 

It follows directly that: 

{6e} - [ G ] { cS A } (III-26) 

Substituting equation IIT-26 in the second term in equation III-15 and 
introducing the notation for the cylindrical element produces the fol- 
lowing expression for the initial stress stiffness matrix, [Ke]. 

[K^ ] - jf {6A} T [G] T [K 1 1 [G]{fiA}Rd0dK (III-27) 

area 

The triple matrix product [ G ] " [K' ] f G j is designated as [HJ and is 
given in Table B-l in Appendix B. The terms in Table B-l are designated 
h^j-j where the i and j denote the row and column, respectively , . in 
which the term is located in the matrix. The overall arrangement of the 
matrix is identical to that shown in Table III-l, and the element corner 
displacement vector { 6 A } has been reorded accordingly. Only non-zero 
terms are given. The matrix is also symmetrical. 

Prior to carrying out the indicated integration, an assumption 
must be made with regard to the character of the stress resultants Np, 
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N 2 and N^2* An assumption, consistent with the assumed linear displace- 
ment functions, is that these stress resultants are constants. With 
this assumption and after carrying out the integration of the right hand 
side of equation III-27, the non-zero terms of the initial stress matrix 
are as given in Table III-3. The arrangement of the terms in the [Ke] 
matrix is the same as indicated in Table III-l. 

Attention is now turned to the third integral in equation III-15. 
Since the displacement functions stated in equations III-21, III-22 and 
III-23 are to be used in the development of the initial dispalceraent 
stiffness matrix, the equation for the initial displacement stiffness 
matrix is obtained in exactly the same manner as was the equation for 
the initial stress matrix, equation III-27. 

Hence : 

[K e 2 ] « ff {6A} T [G] T [K 2 ] [G]{6A}Rd0dX (III-28) 

area 

T A 2 

The triple matrix area product [Gj [K ] [G] is designated as [E] and is 
given in Table B-l in Appendix B. The terms in Table B-2 are designated 
e where the i and j denote the row and column, respectively, in 
which the term is located in the matrix. Again, the overall arrangement 
of the matrix is the same as that shown in Table III-l. The ordering of 
the vector of element corner displacements is changed accordingly. Since 
the matrix is symmetrical, only the diagonal and lower triangular terms 
are given. Also, only nonzero terms are given. 

Before the indicated integration can be carried out, an assump- 
tion must be made about the rotations 3^ and 3 2 - An assumption that 3q 
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and 32 are constants is compatible with the linear displacement functions 
assumed for w. The values of 3-^ and $2 f° r an element are the 
average rotations for that element. The equations used to compute 3^ 
and are given in Chapter IV. With this assumption, the non-zero 

terms of the initial displacement stiffness matrix, obtained from the 
term by the term integration of the right hand side of equation III-28, 


are given in Table III-4; 
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7,2 : 


TABLE III-3 Elements of [K£] 
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TABLE III-3 
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TABLE III-4 ELEMENTS OF [K“] 
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TABLE III-4 (Continued) 
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TABLE III-4 (Continued) 
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TABLE III-4 (Continued) 
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TABLE III-4 (Continued) 
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CHAPTER IV. 

METHOD OF COMPUTATION 

In this chapter, the procedures required to solve for the buckling 
behavior of a complete cylindrical shell structure will be presented. 
Consider the idealization of a complete shell structure to be formed 
entirely of the cylindrical elements developed in Chapter III. In order 
to assess the buckling behavior of the structure, the second variation 
of the strain energy of the structure is formulated in terms of the 
element stiffness, initial stress and initial displacement stiffness 
matrices developed in Chapter III. The concept of combining element 
stiffness matrices to produce a u master stiffness matrix for the com- 
plete structure is well documented (references 27 and 32). 

The second variation of the potential energy for the complete 
structure with the boundary conditions applied has the form; 

5 2 tt = {6i} T Q[K°) + [K 1 ] + [K 2 f){6A} (IV-1) 

where [K°], [K 1 ] and [K 2 ] are the reduced master stiffness, initial 
stress and initial displacement stiffness matrices respectively. 

Applying the stability criterion stated by equation II-2 to 
equation IV-1, the requirement for neutral stability of the complete 
structure is : 

{ 6 A } T Qk°] + [K 1 ] + [K 2 ]]{6A} = 0 (IV-2) 


In the solution of equation IV-2, consideration will be restricted to 
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a linear elastic instability analysis. As used herein, linear elastic 
instability analysis is defined as the calculation of the condition for 
the bifurcation of equilibrium without the need for an iterative deter- . 
mination of the internal loads or deformations. Hence, the magnitude 
of the initial stress stiffness matrix corresponding to the application 
of the critical load is proportional to its magnitude corresponding to 
the application of smaller, but otherwise arbitrary, load. More simply 
stated, 

[K 1 ] I = MK 1 ] ! ; N < N cr 

«cr N 

where the engenvalue, X, is 


A = Ncr 
N 

in which N cr is the load for bifurcation of equilibrium and N is an 
arbitrary initial load. 

In Chapter II, the reason for neglecting the squares of the rotation 
? 

in the [K"] matrix for the beam column was discussed. The same line of 
e 

reasoning, is applicable to the rotations, 3l and 32 » °f the cylindrical 

- 2 , 

shell element. Consequently, the squares of rotations in the [K j matrix 
will be neglected. The initial displacement stiffness matrix for a 
linear elastic stability analysis is 


[K 2 ] 



X [K 2 ] 


N 


The linearized, nontrivial solution of equation IV-2 is 
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det | [K°] + A [[K 1 ] + [K 2 )"] | = 0 (IV-3) 

For computational purposes, the above equation is more conveniently 
stated as 

det | - [I] + [K° ] _1 [[K 1 ] + [K 2 ]] | = 0 (IV-4) 

X — 

where [I] is the identity matrix. 

t 

The major steps involved in the formulation and solution of 
equation IV-4 are: 

1. From basic elemental data, compute the element stiffness 
matrices and construct the master stiffness matrix > [K°] , for the 
structure. 

2. Apply the appropriate boundary conditions to [K°] to form the 
reduced master stiffness matrix, [K°]. 

3. Multiply an arbitrary initial loading by the inverse of [K°] 
to determine the initial displacements, {A}. 

4. Compute the membrane stress resultants, N]_, N 2 and N]_2 for 
each element using 

in) = [s e ] { A e } 

where [S e ] is the "element stress matrix" which is discussed in detail 
and shown in Tables III-7, III-8 and III-9 of reference 25. {A e } are 
the nodal displacements for the element under consideration. 

5. Compute the "average" rotations about the X and 0 axes, $2 


and $1 respectively, from 
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01 = _i [w x - w 2 + - W 3 ] 

and 

3 2 = — + w 9 - w^] + -i [v 1 + v 2 + v 3 + v^] 


6 . Using the results of steps' 4 and 5, construct initial stress 
and initial displacement stiffness matrices for each element and con- 
struct the master initial stress and master initial displacement stiff- 
ness matrices, [K-*-] and [K^], respectively. 

7. Apply the same boundary conditions to [K^] and [K^] as were 
applied to [K°] to determine the reduced master initial stress and 
initial displacement stiffness matrices, [K^] and [K^] , respectively. 

8 . Compute 


[K°] 


-1 



+ 



and compute the eigenvalues of equation IV-4. 

A digital computer program has been written to accomplish the eight 
steps set out above. The program is coded in Fortran IV language and 
has been used on the CDC 6600 digital computer at the NASA Langley 
Research Center. A listing of this computer program is given in 
Appendix C. The reading of input data; the calculation of element 
stiffness, stress, initial stress and initial displacement stiffness 
matrices; the formulation of master and reduced stiffness, initial stress 
and initial displacement stiffness matrices; and the printing of output 
were all coded directly. Library routines (ref. 33) were used for ma- 
trix multiplication, inversion and eigenvalue determination. 
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A concise flow chart of the program is shown in Figure IV-1. The 
general flow of the program follows the previously discussed eight solu- 
tion steps. The subroutine used for the inversion of [K°] uses Jordan's 
method (reference 34) to reduce [K°3 to the identity matrix [I] through 
a succession of elementary transformations. -When these transformations 
are applied simultaneously to [I] and the load vector, the results are 
[K 0 ]-^ an( j displacement vector. The subroutine REIG of reference 33 
finds the eignevalues of a real, square matrix. The original matrix 
which, in this case, is [K 0 ] - "^^] 4- [K^fjis transformed to upper 
Hessenberg form. The eigenvalues are then found using the QR transform 
of J. G. F. Francis (reference 35). 

Because of the vast amount of storage required to solve a problem 
of practical interest, an overlay procedure was used. In the first over- 
lay, the inverted stiffness matrix; the reduced master initial stress 
stiffness matrix; and the reduced master initial displacement stiffness 
matrix are determined. In the second overlay, [K^] and [K^j are added. 
The resulting matrix is premultiplied by [K 0 ] - -*-. The highest eigenvalue 
is then determined for the resulting matrix. Each overlay uses 300,000 
octal storage locations. The computing time, of course, varies with the 


number of degrees of freedom used. A problem havifig about 250 degrees 
of freedom requires about five minutes of computing time. 
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Figure IV. -1 


Elastic instability program flow chart. 
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CHAPTER V. 

APPLICATIONS OF THE COMPUTER PROGRAM 
The procedures described in Chapter II and the finite element 
developed in Chapter III were applied to the stability analysis of sev- 
eral types of structures by means of the computer program outlined in 
Chapter IV. The results of these analyses are delineated in this 
chapter. In addition, investigations pertaining to the influence of the 
initial displacement stiffness matrix and the importance of the non- 
linear terms in that matrix are reported. The types of structures 
considered were: the beam-column, the arch, the flat plate, and the 

curved panel with and without a cutout. The beam-column, flat plate, 
arch and curved panel without a cutout were studied for the purpose of 
establishing the accuracy of the procedure and the finite element. 

Since no information is available on the buckling of curved panels with 
cutouts, the accuracy of those results can only be inferred from the 
accuracy of the solutions obtained for the other types of structures. 

An Euler column with both ends pin-ended was analyzed using the 
finite element developed herein. The exact solution to this problem is 
(reference 12)'. 


P 


cr 


TT 2 EI 


(V-l) 


The column was modeled using plate elements having a width and thickness 
of 1.0 in., thus giving the proper area moment of inertia for the column 
shown in Figure V-l. In this figure, the improvement in the accuracy of 




56 


the solution is shown as a function of the square of the ratio of the 
finite element length (£) to the total length of the column (L) . The 

obtained by H. C. Martin in reference 6 for the same degree of refine- 
ment. 

J 

The present finite element was applied to the solution for the 
buckling of a square, simply supported flat plate. As indicated in 
Figure V-2, the plate was loaded uniaxially with a uniform line load. 

By utilizing symmetry, it was possible to consider only one quadrant of 
the plate. The plate was modeled using square elements of length and 
width, H. 4, 9, 16 and 25 elements per quadrant were used. The influ- 
ence of the mesh size on the solution accuracy is presented in Figure V-2. 
As can be seen in that figure, the accuracy of the solution converges 
rapidly to the closed form solution which is (ref. 12) J 


percent error for an 


>\2 

-] of 0.0625 is in agreement with the solution 


4tt 2 Eh 3 
12(l-v 2 )h 2 


(V-2) 


Also plotted on Figure V-2 arc the results obtained by Kapur and 


Hartz in reference 9. Both the stiffness matrix and initial stress 


stiffness matrix of the Kapur and Hartz finite element were derived using 


the fourth order displacement function for a thin plate in bending which 
was first presented by Melosh (reference 31). The advantage of using 
the higher order displacement function is seen in Figure V-2 to diminish 


as the mesh size decreases. 


have errors of order 



The results .indicate that both solutions 
Since all prebuckling displacements are in 


the plane of the plate, the initial displacement stiffness matrix is 
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identically zero and hence plays no role in the solution of the flat 
plate buckling problem. 

The computer program of Chapter IV was applied to the simply sup- 
ported arch depicted in Figure V-3. The solution to the arch brings 
into play terms containing the element curvature as well as the initial 
displacement stiffness matrix. The exact solution for the buckling of 
a simply supported thin shell arch subjected to a uniform line load is 
given in references 12 and 36 as: 


Qcr “ 




1 


(IX-3) 


For the arch shown in Figure V-3, q cr = 275 lb. /in. The line load 
required to buckle the arch was computed to be 273 lb. /in. using 12 of 
the present finite elements to represent the arch. 

Cylindrical panels of varying curvature were modeled using the 
present finite element and the buckling load was computed. The panels 
considered had equal dimensions in the circumferential and longitudinal 
directions. The panels were simply supported along all four edges and a 

uniform compressive line load was applied in the axial direction. The 
well known cylindrical shell curvature parameter, Z, was varied between 

1.0 and 10.0 by varying the radius of curvature, R.. 

Classical solutions (references 12, 15 and 37) for the buckling of 
curved panels in the curvature range considered predict a single half 
sine wave buckle in both meridional and circumferential directions. 

Since only one quadrant of the panels was modeled, it was necessary to 
establish that the lowest buckling load did indeed correspond to a single 
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half sine wave. This was accomplished by considering two sets of bound- 
ary conditions for the interior edges of the panel quadrant analyzed. 
Both sets of boundary conditions assumed a single half sine wave in the 
circumfei’ential direction but one set assumed symmetry about the midline 
in the longitudional direction and the other set assumed asymmetry. In 
each case, 25 elements were used to model the quadrant. Table V-l gives 
the results of these analyses. The symmetric solutions give the lower 
buckling load in each case, hence, the finite element solution produces 
a buckling mode shape which is compatible with classical solutions. 




z 

SYMMETRICAL SOLUTION 

ASYMMETRICAL SOLUTION 

0 


6.52 

1 

■ 

6.53 

5 

■ Hitkfl 

6.66 

10 

5.08 

7.02 


TABLE V-l BUCKLING COEFFICIENT FOR CURVED PANELS 


As is established by Rehfield in reference 15, boundary conditions 
play an important role in the buckling of curved panels in the curvature 
range under consideration. The classical solution of references 12 and 
37 for curved panels with simply supported edges involve boundary con- 
ditions for the inplane displacements, u and v, which are incompatible 
with the present, analysis. Specifically, during buckling, u is- assumed 
to be zero along the straight edges of the panel and v is assumed to be 
zero along the curved edges. In reference 15, the solution to the 
buckling of a curved panel with boundary conditions which permit inplane 
boundary displacements which are parallel to the simple supports is 
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presented. It is these results which are used as the basis of comparison. 
Figure V-4 shows the convergence of the present analysis to the classical 
solutions of reference 15 as the number of elements is increased. The 
edges of the panels are simply supported. Panels with curvature para- 
meters, Z, of 1, 5 and 10 were investigated. In each case, the answer 
converges to the classical solution rapidly as the number of elements 
used is increased above nine per quadrant. 

Figure V-5 illustrates the geometry of the simply supported curved 
panel with a cutout which was analyzed using the present finite element. 
The circumferential and axial lengths of the panel are equal. The 
cutout is such that its circumferential and axial lengths are also 
equal. Panels having curvature parameters of Z = 0, 1, 5 and 10 were 
considered. Cutout sizes having ratios of cutout length, a, to panel 
length, L, of .1, .25 and .5 were investigated for each panel curvature 
parameter value. In each case 21 elements were used to represent the 
quadrant analyzed. Since the stress distribution and deformations of 
the panel due to the initial load are the parameters involved in the 
initial stress and the initial displacement stiffness matrices, respec- 
tively, i t is instructive to first consider the stresses and deformations 
which result from the application of a unit axial line load to the panel 
shown in Figure V-5- 

Shown in Figure V-6 are stress distributions between the edge of the 
panel and the edge of the cutout along the 0 axis (i.e. , at X = 0) . This 
is the line of maximum stress concentration. The value of the curvature 
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parameter associated with the stress distributions in Figure V-6 is 5. 
However, the plots shown in Figure V-6 are representative of values of 
Z from 1 to 10. The change in the stress distribution as a function of 
curvature, for the range studied, is negligible. As indicated in 
Figure V-6, the influence of the cutout on the stress field is maximum 
at the edge of the cutout and diminishes sharply at points near the 
simply supported edge. The stress distribution shown in Figure V-6 for 

h ci 

the case of y = .1 and .25 are in good agreement with the results given 

by Savin (ref. 40) for an infinate flat plate with a square cutout. In 

reference 41, Savin has shown that the stress distribution is not greatly 

influenced by curvature for the range of Z considered herein. 

Figures V-7, V-8 and V-9 show the deflections, w, normal to the 

middle surface along the x = 0 axis which result from a unit axial line 
load applied to a simply supported curved panel with a cutout. Figure V-7 
presents the results for the smallest cutout considered, a/L = .1. The 
influence of the cutout on the deflections markedly increase as the 
curvature of the panel increases. This same trend is very much in 
evidence in Figures V-8 and V-9 which show the deformations computed for 
the case of cutout sizes of an a/L = .25 and a/L = .5 respectively. 

Thus, while the stress distribution was found to be insensitive to cur- 
vature changes, the displacements normal to the panel middle surface are 
not. In addition, Figures V-7, V-8 and V-9 show the normal deflection at 
the edge of the cutout increases sharply as the size of the cutout 
increases . 

Figures V-10 , V-ll, V-12 and V-13 illustrate the reduction In the 
buckling parameter K x as a function of the cutout size for panels with 
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Figure V. -9 Normal deflection along X=0 due to a uniform axial load 

of -1.0 lb. /in. 
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curvature parameters of 0. (a flat plate), 1., 5. and 10., respectively. 
These plots show the influence of increasing the cutout size is progres- 
sively more drastic as the curvature is increased. This increased sensi- 
tivity to cutouts as the curvature increases is directly attributable to 
the sharp increase in prebuckling deformation which was shown to occur in 
figures V-7, V-8, and V-9 . For example, when the curved panel with Z = 10 
and a/L = 0.5 was analyzed without using [K ], was 4.70 lb /in. How- 
ever, when the problem was resolved, using [K ], was ^*439 lb/in. 

The larger the prebuckling deformations, the greater the influence of the 
initial displacement stiffness matrix. It should be noted that even in 
this case, the prebuckling deformation, w, is more than an order of 
magnitude smaller than the panel thickness. 

The data used to plot figures V-10, V-ll, V-12, and V-13 was cross- 
plotted to produce figure V-14. In this figure, the influence of curvature 
on the buckling strength of panels with various cutout sizes is shown. 

When the cutout Is small, the buckling strength of the panel is seen to 
increase as the curvature increases. This is the same trend exhibited by 
curved panels with no cutout. However, for panels having a larger cutout 
(a/L=0 . 25 and 0.50) this trend is reversed. For the larger cutout sizes, 
increasing the curvature reduces the buckling strength. There are two op- 
posing trends involved here. On the one hand, curvature tends tp stiffen a 
panel, while on the other, curvature increases the magnitude of the pre- 
buckling deformations. 

It is clear from the preceding that the initial deformation stiff- 
ness matrix becomes increasingly important as the hole size and curva- 
ture increase. This raises a question as to the importance of the non- 
linear terms which occur in that matrix. In order to explore this 



Figure V, -10 Buckling of a flat plate with varying cutout size. 
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question, an iterative approach to the problem was taken. This required 
reformulating the problem. First, all of the terms in the [K z ] matrix 
of Chapter III are retained. The [K 2 ] matrix may now be considered the 
sum of two matrices [K 2 *] and [K 2 "]; where [K 2 ’ ] contains only linear 3 
terms and (K 2 ”] contains only squared 3 terms. Equation IV-4 may thus be 
rewritten as 

± [I] + [K° ] -1 [[K 1 ] + [**'] + X ± _! 0 (V-4) 

where 

X-£ is the ith solution to equation V-4 

Equation V-4 was solved iteratively by making successive computer runs 

for the case of a panel having a Z = 10 and a/L = 0-5. This case was 

chosen because it produces large values of 3^ an ^ ^2 an< * ^ ence should 

be the most sensitive to the use of the 3^ terms. The buckling load 

obtained without using the 3 2 term was 0.43939 lb/in. After three 

iterations, the solution had converged to a value of 0.44019 lb/in. This 

2 

is a change of 0.182 percent. Hence, the influence of the 3 terms in 
the initial deformation stiffness matrix is seen to be relatively small 
for the range of parameters considered in this investigation. However, 
for cases involving larger initial deformations an, iterative solution 
would be desirable. 

While no test data is available to directly substantiate the analyt 
ical results obtained for panels with a square cutout, the test data 
obtained by Tennyson (ref. 18) and the conclusions he drew from it are 
in general agreement with the findings reported in this chapter. For 
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example, no reduction in buckling strength is shown for a ratio of cutout 
radius, a, to the cylinder radius, R, of 0.03 but when a/R is increased 
to 0.08, a 40 percent reduction in buckling strength is found to occur. 
Tennyson cites the prebuckling deformation as being a primary factor in 
the reduction of the buckling strength.. 
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CHAPTER VI. 

CONCLUDING REMARKS 

A cylindrical finite element suitable for the linear stability 
analysis of cylindrical shells has been developed. Energy principles 
and variational methods have led to a problem formulation which lends 
itself to physical interpretations of the governing matrices of the 
finite element. By properly grouping the terms which result from taking 
the second variation of the potential energy of the element, it is pos- 
sible to identify three distinct types of matrices. These three matrices 
are : 

1. the conventional stiffness matrix, [K°] 

2. the "initial stress" stiffness matrix, [K^], which is a function 
of the prebuckling stress distribution. 

3. the "initial displacement" stiffness matrix, [K 2 ], which is a 
function of the prebuckling deformations. 

With the assumption of linearity, the buckling problem was stated 
in terms of the classical linear real eigenvalue equation. While the 
stiffness matrix was previously derived, the formulation of the initial 
stress and initial displacement stiffness matrices is orignial. A 
computer program coded in Fortran IV language was developed for use on 
the CDC 6600 series computer. 

The computer program was used to solve several classes of problems 
which have known closed form solutions. Agreement between theoretical 
and computer solutions for the column, the flat plate, the arch and the 
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curved panel are good. The arch solution bears special note since the 
loading in that case is radial. Hence, the applicability of the tech- 
nique to pres sure— loaded structures is assured. A major difficulty 
encountered in the development of the computer program was providing for 
enough degrees of freedom to allow adequate characterization of the pro- 
blem. The computer core storage required by the program is substantial. 
In order to accommodate 36 grid points, representing 216 degrees of 
freedom, overlay programing procedures had to be followed in addition to 
utilizing the entire core storage of 300,000 octal locations. The 
analyses presented in Chapter V indicate that better accuracy could be 
obtained by using more elements. 

The application of. the computer program to the buckling of curved 
panels with cutouts reveals interesting trends. While test data has 
established that, for certain sized cutouts, the buckling strength of a 
cylinder is reduced as the curvature increases; intuition dictates that 
for small cutouts, the stiffening effects of increasing curvature should 
outweigh the detremental effects of a cutout. The analytical results 
for the case of a/L = .1 confirms intuition. 

A number of areas for additional research are apparent as a result 
of this work. By adding beam elements to the existing program, it would 
be possible to evaluate the size of doublers that should be placed 
around the cutout in order to develop higher buckling strength. The 
convergence of the analysis could be improved in several ways. For 
instance, a more complex displacement function could be used in the 
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development of the"initial stress" and "initial displacement" stiffness 
matrices. While this is conceptually straightforward, the work involved 
is formidable. The utilization of an iterative type of solution similar 
to that presented in Chapter V, is another possibility for improving con- 
vergence. Since the first variation of the potential energy, which is the 
basis for the static analysis, actually contains the initial stress and 

i 

initial dispalcement stiffness matrices, the load could be applied incre- 
mentally until bifurcation occurs. These schemes would substantially 
increase computing time. 

Since cylindrical structures with cutouts frequently occur in the 
design of aircraft and space vehicles, it would be highly desirable if 
test programs were initiated to substantiate the analytical findings 
presented herein. 
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CHAPTER X. 

APPENDIX A - THE DEVELOPMENT OF THE CONVENTIONAL , 

THE INITIAL STRESS AND THE INITIAL DISPLACEMENT STIFFNESS MATRICES 
FOR THE BEAM-COLUMN AND ARCH ELEMENTS . 

Development of the Stiffness, Initial Stress and Initial Displace- 
ment Stiffness Matrices for an beam-column and arch are presented in 
this appendix. 

X— 1 The Beam Column 

The equation for the second variation of the potential energy for 
the beam-column element shown in Figure II-l is: 


0 7T = 


f [de 6<j> 6w,J\EA 


10 0 
0 0 0 
0 0 I/A 




0 0 0 
0 N 0 
0 0 0 
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>dX (A-l) 
5w *xx J 


The strains £, <f> and w, xx are related to the displacements u and w 


as follows: 

{if} = [A] {d } 
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{d} is related to u and w through a matrix of differential operators. 


{d } = [D] {g} 


(A-3) 



where 


[D] 
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{g} T = I_u wj 


(A-4) 


Linear and cubic displacement functions for u and w are written as: 


{g} = [b](y) 


where : 
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and { y} is related to the displacements at ends 1 and 2 of the beam- 


column element by 
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Substituting A-3, A-4 , and A-5 into A- 2 


{e} = [A] [D] [ B ] [ r ] {A} 


(A-6) 


Thus 


{6s> = [G]{6A> 


(A- 7) 


where 


[G] = [A][D][B][r] 


Equation A-l may now be conveniently written as 
L L 

6 2 tt = / {6A} T [C] T [K°][G]{6A)dX + / {5A} T [G] T [K 1 ] [G]{6A}dX 

0 0 


+ / {6A} T [G] T [K 2 ] [G] {<5A}dX (A-8) 

0 


where, as in section VI 


[K°] = EA 
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0 <J> o“ 
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; Ik 1 ] “ 

0 N 0 

; [K 2 ] = EA 
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0 0 I/A 
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0 0 0 


< 

Assuming E,A and I to be constants, the first term of equation A-8 


becomes 


{6A} T [K°]{<5A} 




This is the stiffness matrix for a two dimensional beam-column ele 
ment . This matrix agrees with the stiffness matrix for a beam-column 
published in reference 5. 

The second term of equation A-8, assuming N a constant, is 

{so T [K|]ot,) 

where 
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This is the initial stress stiffness matrix for a two dimensional 
beam-column element. This is exactly the matrix obtained in reference 5 
for a beam column initial stress stiffness matrix. 

Assuming <J> to be a constant over the length of the element, the 
third matrix of equation A-8 becomes 

{<5A} T [k|J{«A} 

where 



This is the initial displacement stiffness matrix for a two dimen- 
sional beam-column element. 

X-2 The Arch 

The strain displacement relations for the arch element shown in 
figure A-l may be stated in the form of equation A— 2 by redefining the 
terms on the right hand side of that equation as follows 
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The remaining matrices and matrix operations for the arch are 
exactly the same as those for the beam-column . The resulting conven- 
tional, initial stress and initial displacement stiffness matrices are 
as follows: 
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CHAPTER XI. 

K 

APPENDIX B. - THE PRINCIPAL MATRICES IN THE DEVELOPMENT 
O F THE INITIAL STRESS AND INITIAL DISPALCEMENT MATRICES 
FOR THE CYLINDRICAL ELEMENT 

The explicit statement of the principal matrices involved in deriving 
the Initial Stress and Initial Displacement Stiffness Matrices is given 
in this Appendix. The terminology used herein is consistant with that 
used in Chapter III. 
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-R0(x-a) 

J 


R0-b 

-(R0-b> 

R0 

-R0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-(R0-b) 

R0-b 

-R0 

R0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

(x-a) 

-X 

D 

-(x-a) 

R0-b 

-(R0-b) 

R0 

— R0 

0 

0 

0 

0 

0 

0 

0 

0 

(x-a) 

-X 

X 

“(x-a) 

^-(x-a) (R0-b) 

- f (R6-b) 

xR0 

R 


0 

0 

0 

0 

|(x-a) (R0-b) 

- |(R0-b) 
K. 

xR0 

R 

R0 f v 

- ^(x-a) 

-(x-a) 

X 

-x 

(x-a) 

0 

0 

0 

0 

(x-a) 

R 

_ X 

~ R 

X 

R 

(x-a) 

R 

0 

0 

0 

0 

(x-a) 

2R 

_x_ 

2R 

X 

2R 


i <Re - b > 

- i (R6 - b > 

3R9 

2R 

3R0 
" 2R 

-2 

2 

-2 

2 

— 
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TABLE B-l Elements of [G] T [K 1 ][G] 


ho 9 = No (X-a) 2 (R6-b) 2 
2 R 2 


ho o = - N i o (X-a)(R6-b) 2 _ N (x-a) 2 (R9-b j 
1Z R 2 r 


h 8,2 " “ n 2 


X(X-a) (R0-b) 
R 2 


h q o = Ni 7 (X-a)(R6-b) 2 +N X(X-a) (R9-b) 
r z R 


h 14j2 = n 2 

r\ 2 


h ls o = -Nio R8(X-a)(R6-b) _ N x(X-a)(R9-b) 

R 2 R 


h _ m R0 (X-a) 2 (RQ-b) 

"20 , 2 “ n 2 o 

R 


, R0(X-a)(R0-b) ... ( X -a) 2 (R0-b ) 

" 21 , 2 “ N 12 1 +n 2 ~~ 


h 3,3 = Nx(R6-b) 2 +2N 12 (X-a)(R0-b) +N 2 (x-a) 2 


h 8)3 = N 12 X ( R °-b) 2 +N2 jft -a)(R0-b) 
, R R 


h 9j3 - -Nx (R8-b) 2 -n 12 f(X-a)(R6-b) + X(R6-b) ] -N 2 (X-a) 


„ XR0(R0-b) „ R0 (X-a) 

h 14,3 = - n 12 k " n 2 


R 
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h l5,3 
h 20 , 3 
h 21, 3 
h 8 , 8 = 
h 9,8 = 
h 14,8 
h 15 , 8 
h 20,8 

h 21,8 
h 9j9 = 

h l4,9 ; 
h 15,9 : 


TABLE B-l (Continued) 

= -Ni R8 (R0-b) +N 12 [R0(X-a) +X(R0-bj] +N 2 X(X-a) 

R0 (X-a) (R0-b) R0(X-a)2 

- N 12 — r — +N2 r 

= -N^R0(R0-b) -N 12 [(X-a) 2 - (X-a)(R0-b)J ~N 2 (X-a) 2 

„ X 2 (R0-b) 2 
N 2 ' ^ 

R^ 

X(R0-b) 2 X 2 (RQ -b) 

-Nl2 — “ n 2 r— 

„ X 2 R0(R0-b) 

= -N 2 9 

R" 

„ XR0(R0-b) x2(R0-b) 

= N 12 r +N 2 r 

= N2 XR6 (X-a) (R6-b) 

R 2 

„ XR0 (R0-b) „ X(X-a) (R0-b) 

” “ n 12 ~~ ^ n 2 £ ~ 

N 1 (R0-b) 2 +N 12 (R0-b) (X-a) +N 2 X 2 

= N XR0(RO-b) Ae 

R ^ R 

= -N 1 (R0)(R0-b) -N 12 X[R0+(R0-b)] -N 2 X 2 
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TABLE B— 1 (Concluded) 

,, _ M R0(X-a)(R6-b) M R0X(X-a) 

h 20,9 - ~ N 12 r " N 2 r 


b 21,9 53 N]_R0 (R0-b) +N 2 [XR0+ (X-a)(R0-b)] +N 2 X(X-a) 


u _ M (R9) 2 x2 

h 14 , 14 = N 2 ~ 

R“ 


, „ X(R0) 2 (X-a) 

"20,14 = " N 2 2 


R 


„ X(R0) 2 „ X 2 R6 

h 15,14 = " N 12 — N 2 — 


l 21,14 = n 12 


X(R0) 2 i \t XR0 (X-a) 


+N 2 


R 


h 15 15 = N x (R0) 2 + 2N 12 XRG + N 2 X 2 


(R0) 2 (X-a) R0X(X-a) 


20,15 " 12 . R 


R 


h 21,15 = - N l( R0 ) 2 -N i2 [XR6 + R0 (X-a) ] -N 2 X(X-a) 


h 20, 20 ~ N 2 


(R0) 2 (X-a) 2 


h 21,20 ” “ n 12 


(R0) 2 (X-a) . RO(X-a) 


-N 2 


R 


h 21 21 * N 1^ R8 > 2 +2N 12 R0 ( x " a ) +n 2 < x “ a ) 2 



e ,2,l = 
e 3,l " 
e 8,l = 
e 9 ,1 = 
e l4 ,1 
e 15 ,1 
e 20 ,1 
e 21 ,l 
e 2 , 2 = 


TABLE B-2 Elements of [G] T [K 2 ][G] 

+ ” (X-a) (R0-b ) 2 + ^ (X-a) 2 (R0-b) 

K K 

- (R0-b) [C 2 i(R0-b)+C 51 (X-a>] -(X-a) [C 23 (Re-b)+C 53 (X-a) ] 

- fl 1 - X ( R0 -b ) 2 - — X(X-a) (R0-b) 

R R 

(R0-b) [C 21 (R0-b)+C 51 X] + (X-a) [C 2 3 (R 0 -b) + C 53 X] 

= + XR0 (R0-b) + XR0 (X-a) 

= -(R0-b) [C 2 iR 9 +C 51 X] (X-a) [C 23 R 0 +C 53 X] 

= - R0 (X-a) (R0-b) - -=~ R0 (X-a) 2 
= (R 8 -b) [C 2 iRe+C 31 (X-a) ] + (X-a) [C 2 3 R 0 +C 53 (X-a) ] 

- (X-a) (R0-b) 2 + (X-a) 2 (R0-b) 

R R 

1 C 55 

+ i. (X-a)(R0-b)[C 35 (R0-b) + C 45 (X-a) + (X-a)(R0-b)] 

R K 


e 3>2 = -(R0-b) [C 23 (R0-b) + C 53 (X-a)] - (X-a) [C 24 (R0-b)+C 5 A (X-a) j 
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TABLE B-2 (Continued) 

C o 

e 7,2 = “ — ~(X-a) (R0-b) 2 -C 53 | (X-a)(R0-b) 

C 53 X 7 c 54 x 

e 8>2 = - “Y" < Re " b > " ”5“ (X-a)(R0-b) 

- i (X-a) (R0-b) [C 35 (R0-b) + C 45 X 4- C 55 j (R0-b) ] 

R 

e 9,2 = (R 6 -b) {C 23 (R0-b) + C 53 X] + (X-a) [C 2 4 (R 0 -b) + C 54 X 3 
+ | (X-a)(R0-b) [C 25 (R0-b) - C 45 £ (R0-b) + C 55 X] 
e 13,2 " C 51 Y ( X-a )( R 6 - b) + C 53 (X-a) (R0-b) 

e 14,2 “ C 53 ™ (R9-b) + C 54 Y (X-a)+ i (X-a) (R0-b) [C 35 R0+C 45 X+C55™] 
e 15,2 = -( R 0“b) [C 23 R 0 +C 53 X 3 - (X-a) [C 24 R 0 +C 54 X] 

+ ^ (X-a)(R9-b)^C 25 R +C 45 - C 55 X] 

e l9 1 2 = ” c 51 y < x “ a ^ R8 - b > " "Y ( x -a) 2 ( R Q-b) 
e 20,2 = " C 53 y ( Re ” b ) -C 54 ( x " a > 2 

5 K 

- |(X-a)(R0-b) [C 35 R0+C 45 (X-a) + C 55 M (X-a)] 
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TABLE B-2 (Continued) 

e 2 i 2 = (Re-b)[C23R0+C 5 3(X-a)] + (X-a) [C 24 Re+C54(X-a) ] 

+ i (X-a)(R9-b)[C25R0-C 4 5 ^ (X-a) + C 55 (X-a)l 

2 

e 3 3 = C 22 (R0-b ) 2 -C 42 ^ (X-a) (R0-b ) 2 +C 52 (X-a) (R0-b) -C 2 4 ^ (X-a)(R0-b) 

C c o 

— (X-a) 2 (R0-b) +C 25 (X-a) (R0-b) (X-a) 2 (R 6 -b) +C 55 (X-a) 

R “ 

e 7 3 « (R0-b) [C 2 i(R 0 “b) + c 5 l( X-a )l + [C 2 3 (R 0 -b) +C 53 (X-a) ] 

e 8 3 3 023^®“ b) 2 +C 53 (X-a)(R0-b) + C 24 X(R0-b) + C 34 X(X-a) 

+ C 25 f (R0-b) 2 - C 45 (X-a) (R0-b ) 2 + C 55 |(X-a)(R0-b) 

r 2 R 

0 y 0 ^24- 2 

e 9,3 = -C 2 2 (R^b)^ + C 42 | (R0“b) z - C 52 X(Re-b) + — (X-a)(R0-b) 

+c 54 f (X-a) (R 8 -b) -C 23 (X-a) (R0-b) +C 4 5 § (X-a) (R0-b) -C 5 5 X(X-a) 

e 13 3 = - c 21 R0 ^ Re “ b ^ “C 51 R0(X-a) -C 23 X(R0-b) -C 53 X(X-a) 

e 14 3 = “ c 23 R0 ( R0 " b ) -C 53 R0(X-a) -C 24 X(R0-b) -C 54 X(X-a) 

+ X R 0 [-C 25 (R 6 -b) +C 45 I (X-a) (R0-b) -C 55 (X-a)] 

R R 
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TABLE B-2 (Continued) 


e 15,3 = C 22 R0(Re-t>) -C 42 SE (R0-b) +e 52 X(R 8 -b) - C 2 4 “ (X-a)(R0-b) 

-C 5A | (X-a) (R0-b) +C 2 5 R 6 (X-a) -C 45 (X-a) +C 55 X(X-a) 

e 19 } 3 = C ? 1 RO(R0-b) +C 51 R0(X-a) +(X-a) [C 23 (R0-b) +C 53 (X-a) ] 
e 20,3 = c 23 r0 ^ r +C 5 3 R 0 (X-a) +C 24 (X-a) (R0-b)+ C 5 4 (X-a ) 2 

+C 25 M (X-a) (R0-b) -C 45 ~ (X-a ) 2 (R0-b) +C 55 M (X-a) 2 
e 2l,3 = -C 22 R0(R6-b) +C 42 M (X-a)(R0-b) -C 52 (X-a) (R0-b) +C 24 (X-a)(R0-b) 

C <> 4 9 R0 /s 9 

+ — (X-a) z (R0-b) -C 25 R 0 (X-a) +C 45 -pr (X-a ) 2 -C 55 (X-a) z 

C 2 

eg 7 = -# x ( R0 " b ) 2 + c 53 ^ (R 0 - b ) 

eg } 7 = -(R 8 -b) [C 21 (R0-b) +C 51 X] -X[C 23 (R0-b) +C 53 X] 

X0R X 2 0R 

e 14 , 7 = " c 51 R (R0 _ b) -C 5 3 R 

e 15,7 = (R0~b) [C 2 xR0+C3iX] + X[C 23 R0+C5 3 X] 

fiR R0X 

e 20,7 = C 51 T (X-a)(R0-b) +C 53 (X-a) 

e 2l,7 = ~ (R 6 -b) [C 21 R0+C 51 (X-a) J -X[C 23 Re+C 53 (X-a)] 
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TAM/E B-2 (Continued) 


eg, 8 = (-'53 f (HO-b) 2 +C54 (RO-b) +C35 | (R8-b) 2 


2 v2 

+C45 (RO-b) +C’5«- > ~ ; y (RO-b) 

R- 


e 9j8 ■= j ( KO-L) ) 2 -C r> jX(KO-b) -C24 X ( RO-b) -C54X 2 -C 2 y ^ (RO-b) 2 


+C 4 , JL (RO-b) ^ -C53 ---- (RO-b) 

2 R 

R 


2 


X 2 


e 13 8 = ~ c 51 X -°- (KO-b) “O53 --- (RO-b) 
R R 


„ XRO ,, jn , , r , X 2 RO XRO ,_ n . v 

G 14 , 8 = “ c 53 — B - ( K °-b) -C54 ~ ‘ r“ ( R °- 0 ) 


X 2 1 X 2 R0 

-C45 -V (RO-b) - C 55 ----- (RO-b) 


e 15,8 = (’2 3R (RO-b) +C53X ( RO-b ) +C24ROX +C54X 2 +C25 ----- (RO-b) 

-C 45 (KO-b) +<>>5 X O(0-b) 

e J.9,8 = O5J -yp (RO-b) +C53 ~ (X-a) (RO-b) 

e 20,8 = O53 (X-a) (RO-b) + C 54 -jp- (X-a) +C35 ~~ (RO-b) 

+C45 ~ (X-a) (RO-b) +C55 XR y- (X-a) (RO-b) 

K R z 
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TABLE B-2 (Continued) 

e 2l,8 = -C23RQ(R0-b) -C 53 (X-a) (R0-b) -C 2 4 xr9 -C 54 X(X-a) 

-C 2 5 ■— <R0-b) +C 45 ™ (X-a) (R9-b) -C 55 £ (X-a) (RQ-b) 

R R 

e 9,9 “ C 22 (R0-b ) 2 “C 42 | (R0-b ) 2 +C 52 X(R0-b) -C 2 4 f (R0-b ) 2 

2 • 2 
-C 34 \ (R 6 -b) +C 25 X(R0-b) -C 45 — (R0-b) +C 55 X 2 . 

R R 

e 13,9 = C 2 lR0(R0-b) +C 51 R0X +C 2 3 X(R 6 -b) +C 53 X 2 
e 14,9 = C 2 3 R 0 (R 6 -b) +C 33 XR 0 +C 24 X(R0-b) +C 54 X 2 +C 25 (R0-b) 

-C 45 (RQ-b) + C 55 ^ 

Y A 0 

e 15,9 = ” c 22 r ( R 0-b) +C 42 ( R 0“b) -C 52 X(R0-b) +C 24 ( R 0-b) 

x z R 0 

+c 54 "r* ( R 0“b) -C 25 XR 6 +C 45 -C 55 X 2 

e l9 9 * -C 2 i r 0 (R0-b) -C 51 X0R - (X-a) [C 2 3 (R 0 -b) +C 53 X] 

e 20,9 = -9 23 r ® ( R 9_ b) -C 53 XR 8 -C 24 (X-a) (RQ-b) -C 54 X(X-a) 

-C 25 M (X-a) (R0-b) +C 45 (X-a) (R0-b) -C 55 ~ (X-a) 

R 
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TABLE B-2 (Continued) 

ufl XR 8 

e 2l,9 = C 2 2 R 8 (R 6 -b) -C 42 X (X-a)(R0-b) +C 52 (X-a)(R 6 -b) -C 2 4 


-C 54 ~ (X-a)(R0-b) +C25XRO “C45 (X-a) +C 5 5X(X-a) 


XR 2 6 2 . X 2 0R 

R 


e 14 , 13 *= C 51 ' R +c 53 R 


e 15, 13 * -C 2 iR 2 0 2 -C 51 XR 6 -C 23 XR0 -C 53 X‘ 


R 2 0 2 „ X0R /v x 

e 20,13 C 51 R ( x_a ) ~ c 53 nT - 


e 21 13 = c 2lR 2 ^ 2 +C 51 (X-a)R0 +C 2 3XR0 +C 53 X(X-a) 


XR z 0 


2 0 2 „ X 2 R0 XR 2 0 2 _ X 2 R0 X 2 R 2 0 2 


e 14 , 14 = c 5 3 — r — +c 54 +g 35 r"~ + c 45 r*" fC 55 


R 


2n ? „ „ XR 2 0 2 

e 15,14 = -C 23 R^ -C 33 XR 0 -C 94 XR 0 -C54X 2 -C 25 — r— 


r XR 2 0 2 XR0 v 

e 19 , 14 = " c 51 c 53 ~r~ ( x " a ) 


„ R 2 0 2 ^ XR0 , , N „ XR 2 0 2 „ XR0 . . 

G 20 , 14 = -C53 ~r (X-a) -C54 ~ — (X-a) -C35 — r -C45 -y- (X “ a) 


-C55 M?i - 2 (X-a) 
R 2 


(R0-b) 



107 


e 2l ,14 


e 15 , 15 


e 19 , 15 
e 20 , 15 


e 21 ,15 


e 20 , 19 
e 21 ,19 
e 20 , 20 


TABLE B-2 (Continued) 


c 23 R 202 


+C 5 3R0(X-a) 


+C24XR0 +C 34 X(X-a) 


+ C25 


XR 2 6 2 

R 


” c 45 


XR 2 0 2 


R 


2 


(X-a) 


+c 55 


XR0 

R 


(X-a) 


^ 22 r2 ® 2 “^42 


XR 2 6 2 

R 


XR 2 6 2 X 2 R0 

+C 52 XR0 -C 24 ^ — -C 54 •— g— +C 25 XR 0 


-C 45 ^ + C55X= 


C 2 iR 2 0 2 +C 51 R0X + (X-a)fC 2 3 R 0 +C 53 X] 


>2 a 2 


R 2 © 2 


c 23 r “9 +C 53 XR0 +C 24 RO (X-a) +C 5 4 X(X-a) +C 25 — — (X-a) 


-C 45 ^ (X-a) +C 55 -T (X-a) 


-C 22 R 2 9 2 + C 42 (X-a) -C 52 R0 (X-a) + C 24 — -- 

K . R 


XR0 XR0 

+ C 54 — (X-a) -C 2 5 XR0 +C 45 (X-a) - C 53 X(X-a) 


r r2q2 (v xr R0 /v n2 

C 51 ~ jp~ (X-a) +C 53 -g- (X-a) 


-C 2 iR 2 0 2 -C 51 R0(X-a) - (X-a) [C 23 R9 +C 5 3 (X-a)] 


= C 


r 20 ' 


53 


, R 0 ^ 


R 2 © 2 


R0 


(X-a) +C 54 R (X-a) +C 35 " - • (X-a) + 045 —— (X-a) 


+c 55 (X-a ) 2 
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TABLE B-2 (Concluded) 

2 r2q2 

e 21 , 20 = ~^23 r2 ® 2 “C 53 R 8 (X-a) -C 24 R 8 (X-a) -C 54 (X-a) -C 25 — (X-a) 

R 2 8 2 .2 „ R 8 , v *2 

+C 45 ~2~ ( x “ a ) “ c 55 < x ~ a > 

R 

o o r 2 o 2 d 2 q 2 R 0 . .2 

e 21,21 = c 22 r2 ® 2 " c 42 ( x ’ a > +c 52 R0 ( X-a ) _c 24 (X-a) _C 54~Y ( X-a ) 

+C 2 5 R 0 (X-a) -C 45 ~ (X-a ) 2 +C 55 (X-a ) 2 
R 


where : 


C 21 - Pi; C 51 = V S2 ; c 22 = B l 2 + (1 2 V> 622 


C 32 = ^ $ 2! c 42 - «Bi: C 52 - h®2 


C 54 - »i: c 55 “ e 2 2 + Pi 2 


2 
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CHAPTER XII 


APPENDIX C - COMPUTER PROGRAM LISTING 
The computer program used to obtain the results presented in 
Chapter V is presented in detail in this Appendix. This program is 
called STABL. The input data required are as- follows; 

NC number of cases to be run 

number of finite elements 
N <f)DE number of nodes 


NE 


AZ(MN) 

BZ(MN) 

CZ(MN) 

EZ(MN) 

TZ(MN) 

XMUZ(MN) 


Nl(MN), N2 (MN)1 


N 3(MN) , N4(MN)j 


JR(I) 

f4>rc(i,1) 


number of degrees of freedom to be constrained 
length of the finite element, MN, in the X-direction 
length of the finite element, MN, in the 6-direction 
curvature of the finite element, MN 
modulus of elasticity of finite element, MN 
thickness of finite element, MN 
Poisson's ratio for finite element, MN 
the four mode points of finite element 
MN, read counter clockwise, with N1 and N2 estab- 
lishing the element X-axis 

a list of the degrees of freedom to be restrained 
the vector of applied forces 


A listing of the program follows. Comment cards are included in 


the listing to provide clarifications of program functioning and 
t erminology . 
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OVERLAY (LINK *0,0) ♦ 

PROGRAM STABL ( INPUT, OUTPUT * TAPE2 . TAPE3 ) 

COMMON/ZZ 1 /AZ ( 30 1 *BZ < 30 ) « CZ ( 30 > * DMZ ( 30 > « DBZ < 30 > * XMUZ ( 30 > 

1 /ZZ2/XK0(24*24 > * 

2 /ZZ3/N1 (50 ) *N2 (50)«N3(50) *N4 ( 50 ) 

4 /ZZ5/JRU50) 

5 /ZZ6/FORC < 21 6 * 1 > 

6 /ZZ7/X (216) 

7 /ZZ0/S (3*24 ) • 

8 /ZZ9/STRSR ( 3 > * 

9 /ZZ10/XK1 (24*24 ) . 

9 /ZZ 1 1 /XE ( 24 ) 

9 /ZZ1 2/BTA1 «BTA2 

9 /ZZ I 3/XK2 ( 24 t 24 > 

9 /ZZ 1 4 /SO 

L I NK =: 4 LL 1 NK • 

read 1000 »NC * SO 

1000 FORMAT <I3*E12*4) 

N0C = 0 

1001 CALL OVERLAY (LINK* ! *0*0 ) * 

CALL OVERLAY (LINK «2i 0*0) ♦ 

N0C=N0C+1 

IF(NOC*LT*NC> GO TO 1001 
STOP « 

END* 

OVERLAY (L INK* 1*0)* 

COMMON/ZZ 1 /AZ ( 30 ) *BZ ( 30 ) * CZ ( 30 > * DMZ (301* DBZ ( 30 ) « XMUZ (30 ) 

1 /ZZ2/XK0 (24*24 ) • 

2 /ZZ3/N1 (50 ) ,N2 (50 ) *N3 (50 > *N4 (50 ) 

4 /ZZ5/JR < 1 50 > 

5 /ZZ6/F0RC (216* ! ) 

6 /ZZ7/X (216) 

7 /ZZ8/S ( 3,24 ) « 

8 /ZZ9/STRSR ( 3 ) ♦ 

9 /ZZ10/XK1 (24*24 1 ♦ 

9 /ZZ 1 1 /XE ( 24 ) 

9 /ZZ12/BTA1 oBTA? 

9 /ZZ1 3/XK2 (24*24 ) 

9 /ZZ I 4/SO 

D I MENS ION IP I VOT (21 6 ) * XMK (216*216)* INDEX! 21 6 * 2 ) * EZ ( 30 ) « TZ ( 30 ) ♦ 

C A IS LEN* IN X DIR* B IS LEN. IN THETA Dl R* C IS CURVATURE * XMU IS 

C POISSONS RATIO* XKO [S THE ELEM STIFF# MATRIX# 

C no=number OF ELEMENTS *' 

C NODE^NUMBER of nodes « 



o u 


111 - 


c IF S6 = 0. BETA SQRD TERMS ARE IGNORED* IF S0=1* BETA SQRD TRMS USEC 

C NE^NuMBER OF DOF TO BE RESTRAINED ♦ 

1001 READ 1 01 *NO*NO0E*NE « 

101 FORMAT (313) * 

J00F=6*N0DE * 

C INITIALIZE MASTER ST IFFNESS MATR I X t XMK * TO ZERO 

DO 51 I 1 = 1 *JDOF • 

DO 51 JJ = 1 • JDOF « 

XMK ( I I ♦ JJ ) = 0 • 0 * 

51 CONTINUE * 

DO 1 MN=! « NO * 

Read 1 05 * AZ (Mn > *BZ (Mn ) * CZ (MN ) «EZ (MN) ,TZ (MN) * XMUZ CMN ) 

105 FORMAT (6E12*4) • * 

DMZ ( MN ) =EZ ( MN )*TZ (MN ) / ( 1 * -XMUZ ( MN ) **2 ) 
DRZ(MN)=DMZ(MN)*TZ(Mn)**2/12. 

C READ THE NODE NOS IN COUNTER CLOCKWISE DIR. FOR ELEM. MN 

READ 109. Nt (MN) «N2 (Mn>*N3(MN> *N4(MN) 

109 FORMAT (413) . 

PRINT 53. MN 

53 FORMAT (/40X .*ELEMENj NUMBER* *13) 

PRINT 54 » AZ (MN) *BZ(Mn) *CZ ( MN ) * EZ (MN) *TZ (MN) *N1 (MN) * N2 (MN ) *N3 ( MN ) • 

1 N4 < MN ) 

54 FORMAT </lX**A=*«E12#4**B=*«E12*4«*C=**E12*4.*E=**E12*4.*T=*.E12.4 
1 ,*N1 =*. I 3 «*N2 = * • I 3.*N3 = * .13* *N4 = *. I 3 ) 

CALL ELEMKO (MN) * 

CALL ADDUP (MN*XMK.X<0) 

1 CONTtNUE « 

C READ IN THE D.O.F. TO BE EL I MINATED* JR( NE ) 

READ 121 * ( JR(1 ) , I =1 *NE) * 

121 FORMAT (2513) * 

PRINT 55 

55 FORMAT </40X.*DEGREES OF FREEDOM TO BE ELIMINATED*) 

PRINT 1 21 • { JR ( I ) * I = 1 ,NE ) 

ELIMINATE NE DEGREES OF FREEDOM FROM XMK(I.J) TO OBTANTAIN THE 
REDUCED MASTER STIFF MATRIX WHICH WILL STILL BE CALLED XMK BUT IS 
C OF ORDER JDOF-NE , 

CALL WASH ( JDOF* NE* XMK) 

JOOFR= JDOF—NE , 

C ZERO OUT APPLIED FORCE JECTOR « 

DO 3 1=1*216 

FORC ( I * 1 ) =0.0 , 

3 CONTINUE « 

C READ IN APPLIED FORCfS * 

READ 2* (FORC ( I . I ) *1=1 .JDOF ) 
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2 FORMAT (6E12.4) 

C REDUCE FORCE VECTOR TO CORRESPOND TO RED* MASTER STIFF* MATRIX 

CALL REDFORC (JDOF*NF>* 

C PRINT REDUCED FORCE VECTOR* STILL CALLING IT FORC 

PRINT 133 ♦ 

133 FORMAT (/40X **FORCES APPLIED TO UNRESTRAINED D.O.F*/) 

PRINT 1 37* C I *FORC (1*1 1*1-1 « JDOFR > 

137 FORMAT ( 4X* *FORC (* * | 3** > = * *El2.4 ) 

C SOLVF 'FOR DISPLACEMENTS • 

CALL MATINV (XMK, JDOFR*FORC*l *DETERM*1PIV0T*INDEX.216.ISCALE) 
XMK- INVERSE IS NOW STORED IN XMK 

NAME DISPLACEMENTS. AT THIS POINT THEY ARE STORED IN FORC 
DO 5 1 = 1 * JDOFR « 

X< I )=FORC (1*11 
5 CONTINUE ♦ 

REORDER DISPLACEMENTS TO AGREE WITH NODE NUMBERING 
CALL EXPDEF <JDOF*NE>* 

PRINT OUT DISPLACEMENTS ♦ 

PRINT 141 • 

PRINT 145* < I *X< ! 1 * 1 = 1 * JDOF 1 ♦ 

REWIND 2* 

WRITE TAPE 2 • ( ( XMK < I * J 1 ♦ I = 1 * 1 SO ) * J = 1 « 1 80 ) «NO» JDOF »NE * JDOFR* 

141 FORMAT (/50X**DtSPLACEMENTS*/> 

145 FORMAT ( 4X**X(*t 13* * >=+*E12*4 ) * 

DEVELOP THE STRESS MATRIX FOR EACH ELEMENT AND COMPUTE THE* 
INPLANE STRESS RESULTANTS AT THE CENTER OF THE ELEMENT* 

DO 1 3 I =1 « 21 6 
DO 13 J=1 *21 6 

I 3 XMK ( I * J ) = 0 • 

DO 11 MN= 1 » NO * 

CALL STRESS (MNl* 

COMPUTE THE ELEMENT PRESTRESS MATRIX*XK1 
CALL ELEMK1 (MN) * 

COMBINE ELEMENT XK1 MATRACIES IN TO A MASTER PRESTRElSS MATR I X 
CALL AODUP <MN*XMK*XKll 

COMPUTE ELEMENT PREBUCKLING DEFORMATION MATRIX 
CALL ELEMK2 (MN) 

COMBINE ELEMENT X<2 MATRACIES INTO A MASTER K2 MATRIX 
CALL ADDUP(MN*XMK.XK2) • 

II CONTINUE* 

THE MASTER <1 MATRIX IS IN XMK 
REDUCE XMK l TO PROPER ORDER 
CALL WASH ( JOOF*NE*XMK > 

REWIND 3 
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WRITE TAPE 3» ( ( XMK (I»«J>.!®1* JDOFR ) * 1 . JDOFP ) 

RETURN* 

ENO. 

SUBROUTINE ELEMKO <MN) * 

COMMON/ZZ 1 /AZ ( 30 > ,BZ { 30 > . CZ ( 30 > . DMZ ( 30 ) . DBZ 130). XMUZ ( 30 ) 
t /ZZ2/XK0(24.24 > * 

2 /ZZ3/N1 ( 50 ) *N2 ( 50 ) . N3 ( 50 ) *N4 < 50 ) 

A=AZ(MN> 

b=bz (MN ) 

C=CZ(MN) 

DM=DMZ(MN) 

OB=OnZ(MN> 

XMU=XMUZ(MN) 

AB= A/B • 

BA-B/A « 

ATB=A*S . 

AB3“ A/B**3 • 

BA3=R/A**3 « ' , 

XKO < 1 * l )=XK0 <7.7)“XK0 ( 13* 13>“XK0 ( 19.19) “DM# < BA/3.+ AB* < 1 • — XMU 1/6 • ) 
XKO (2. 1 )“X<0 < 14. I3)=XK0C20.7) = XK0( 19.8>=DM#( 1 .+XMU)/8. * 

XKO <3.1 ) “XKO (9.1 )“XKO Cl 9.15) =XKO ( 21*1 9 ) “-7.*XMU*B*DM*C/40 • 

XK 0(4*1 ) =XKO C 1 0 *7 )-XKO <16.13 )*XKO (22 * 19 ) “XKO O * 1 >*A/6. 

XKO (5. 1 >=XKO (11*1 ) “XKO (17.13 )“XKO (23.13) =-B**2*XMU*DM*C/4 0 • 

XKO ( 6. 1 ) “XKO (12.7)= XK O ( 2 A . 13) = XK 6(19.18)“ XK 0(5.1 )*A /6 • 

XKO (7.1 )=XK0(19.13)=-DM*(BA/3.-( 1 .-XMU)*AB/1 2. > 

XKO (8.1 )=XKO< 19.2)=XKO(14.7)=XKO(20.13)=-DM*( 1 *-3.#XMU >/8 » 

XK 0 ( 1 0 . 1 ) -XKO (22. 13) = XK0( 19. 1 6 ) “XKO < 7. 4 ) “-XKO (4.1 ) 

XKO { 1 2. I )s:XK0 (7.6 ) “XKO ( 1 8* 1 3 ) “XKO ( 24 • 1 9 ) *-XKO < 6 » I ) 

XKO < 1 3. 1 )“XKO ( l9.7)“-XK0 (1*1 >/2» 

XKO (14.1 ) “XKO (13.2) “XKO (8*7) -XKO (20*19) =— XKO (2.1 ) 

XK 0(15.1) =XKO (21.1 ) “X* 0(19*3)= XK 0(19.9) “-3.*XMU*B*DM#C/40 • 

XKO < i 6* 1 ) “XKO (1 3.4 >“X<0( 22.7 )“XKO( 1 9 * 1 0 ) “ATB*XMU*DM*C/80 • 

XK 0(17*1 ) =XKO ( 23 . 1 )*XK0<13»1 1 >=XK0<13*5)“ B**2*XMU*DM*C/60 . 

XKO < 18. 1 )=XK0(l9.6) = XK0(24.7)=XK0(l3.12)=XK0(5.i )*A/9. 

XKO ( 1 9. 1 )=XK0 ( 1 3.7)“OM* (BA/6.- ( 1 .-XMU)*AB/6. ) 

XKO (20.1 )=XKO<7.2)=XKO(13.8)=XKO( 1 9 « 1 4 ) =-XKO ( 8 . 1 ) 

XKO (22 i 1 ) =XKO (19.4) =XK0 (16*7) =XKO (13*10) XKO (16*1 ) 

XK0(24. 1 )=XKO ( I 3.6)=XK0( 18.7 )=>XKO ( 19. 12 ) = -XKO ( 18* 1) 

XK 0(2*2 ) = XK 0(8.8) =XKO < 14*14) “X<0 ( 20.20 )= DM* ( ( AB/3. )* ( 1 .+D8* (C**2 ) / 
1 DM )+ ( ( 1 •— XMU )*BA/6. )*< 1 • *4 • *DS* ( C**2 ) /DM > ) 

XKO (3*2 ) =XK0 (9.8 ) =-7.*A*DM*C/40.+DB* (2.-3.*XMU )*C/ ( 2.* A ) 

XKO (4*2 ) =XKO ( 22 *2)=XK 0(14.10) “XKO (16.14) s-«A**H#DM*C/40.-XMU*DB*C/ 

12. 

XKO (5 .2 ) “XKO ( 1 1 .8 )“XKO (17.14)*XK0(23.20 >“-7.*ATB*DM*C/Z40 • -D0*C* ( 
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1 <2*-XMU)*B/< 12. *A >+7«*AB/20. 1 

XKO (6*2 )=XKO ( 1 8*®-) = XKO <14*12)=XKO<24,20 ) =-A*ATB*DM*C/240 .-DB*C* < 

1 XMU*8/12.+A*AB/20. > 

XK 0(8*2 >=XKO (20* 14 } =DM* ( AB/6 • * < l .+DB*C**2/DM >~ ( 1 *— XMU )*BA/6o* < 1 e + 

1 4.*DB*C**2/DM> j 

XKO <9*2 >=XKO<0,3>=~3.*A*DM*C/4O.-<2.-3.*XMu)*DB*C/'<2.*A > 

XK 0 ( i 0 * 2 ) =XK0 < 16*2)=XK0 < 1 4 ♦ 4 > =XK 0 < 22 . 1 4 ) = DM*A**2*C/60 . 

XKO < l l ♦ 2 ) -XKO (B* 5 ) = XKO <23,14>=XK0<20*17 ) = ~ATB*DM#C /SO o +DB*C* < <2*- 
1 XMU 1 *B A / 1 2 ■ — 3 • *■ AB/20 . > 

XKO < 12, 2 )=XKO ( 14 ,6) =X*0 <24*8 > = XKO <20*18 >=A*ATB*DM* C/360 .+A*AB*DB 
l*C/30* 

XKO < 14*2 >=XK0.< 20,8 1=-XK0< 2*2 >/2» 

XKO< 15*2 >=XKO <21 .8)=-3.*A*DM*C/4O.+<2.-XMU)*0B*C/ <2.*A ) 

XKO < i 7*2 )=XKO < 14*5 > = XKO (23*8 > = XKO <20* 1 1 >=-XKO< 1 1 *2 1 

XKO <18,2 > = XKO <8,6 ) = XKO <20* 12 > = XKO <24 . 1 4 > =-XKO ( 12*2 1 t 

XKO<?0, 2 ) -XKO < 14 ,B)=-DM#< AB/3.* < 1 • +DB*C**2/DM } -B A * < 1 .-XMU>/J 2«* < l * 
1+4.*0B*C**2/DM> ) 

XKO <21 , 2 ) =XKO < l5.8)=-7#*A*DM*C/40.-<2«-XMU)*DB*C/(2**A > 

XKO <23,2 )=XKO <20*5 J=XK0 < 17*8>=XK 0(14*11 >=-XK0<5*2) 

XKO <24, 2 >=XKO <20*6 ) -XKO < 1 2 , 8 ) = XK 0 U 8 * 1 4 )=-XK0<6*2 ) 

XKO (3,3 >=XK0(9,9)=XK0< 1 5 , 1 5 > =XKO ( 2 1 *21 ) - 1 56 ./35.*DB* <B/A**3+A/ 
l9*'#3) + 72./25. +DB/ ATB+ 1 69 •* ATB*DM*C**2/ 1225* 

XKO <4* 3 ) -XKO <22*21 1= A* < DB* ( 78 « /35 • *B/A**3+22 ./35* *A/B*#3+l ./ATB*i 
16./25* + 6**XMU/5. ) > + 1 43 .*ATB*DM*C**2/ < 6# * 1 225 • ) > 

XK 0(5*3)= XK 0 < 1 1 *9> = e*(0B*<22.*B/ < 35 .*A**3 > +78 .*A/ <35. *8**3 ) + l «/ 

1 ATB*<6./25.+6.*XMU/5, ) ) + l 43* *ATB*DM*C**2/ <6* * 1 225. } 1 
XKO (6* 3 >=XKO< 18* 15)= A*B*(DB*< 1 1 .*B/ < 35* *A**3 > + 11, /35.*A/B#*3+ 

1 1 •/ATB»< ,02+.2*XMU> ) + l 21 .*ATB*DM*C**2/<36.*1225. > > 

XKO <7.3 >=XKO <9,7>=XK0 <15. 1 3) =XK0<21 ♦ 1 3 1 = -XKO < 3 , 1 > 

XKO <9* 3 >=XKO (21 . 15 ) =-DB*< 1 56 */35 . *8/ A**3-54 * /35 • *A /B**3+72 o / < 25 
1 ATB) > + l 17.*ATB*DM*C**2/<2.*1 225. ) 

XKO < i 0 • 3 ) = XKO (21*16)=: A* < DB* <78./35**B/A*+3-l 3./3S. *A/B**3 +6./< 

1 25.*ATB >1-1 69,*ATB*DM*C**2/< 12**1 225. > ) 

XKO < 1 1 . 3 ) =XKO <9.5 )= B*<DB*<-22./35.*B/A**3+27./35.*A/B**3-l ./ATB* 

1 <6./25.+6.*XMU/5. ) >+33.*ATB*DM*C**2/ <4. *1 225. ) > 

XKO( 12. 3)=XK0 <24* 15>= A*B*(DB*( 1 1 ./35,*B/A**3- 1 3. /70. * A/B**3+l ./ 

1 ATB*< • 02+XMU*. 1 ) )— 1 43.*ATB*DM*C**2/ <72**1225. ) > 

XKO ( 1 3 * 3 ) = XKO ( 15,7) =XKO <21 .7 >=XKO < 1 3.9)=-XK0 < 15*1) 

XKO < 14* 3 ) =XKO (20*9 ) =-XKO < 1 5* 2} 

XKO ( 15* 3)=XK0 <21 .9)=-54./35.*DB*<B/A**3+A/B**3 )+72 . /25.*DB/ATB+ 

1 81 .*ATB*DM*C**2/ <4.*l 225. > 

XKO < 16*3 )=XK0 <21*10)= A* (D0* ( 27 • /35 • *B/A**3+ I 3 . /35 • *A/B** 3-6 • / 

1 <25.*ATBn-39.*ATB*DM*C**2/<8.*1225. > > 

XKO< 17.3>=XK0<23*9>= 8* < DB* < 1 3* /35 .*B/A**3+27 . /35 * * A/B**3-6® / 
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1 (25.* ATB) )-39.*ATB*0M*C**2/(8.*1225. ) ) 

XKO (18*3) =XKO ( I 5 * 6 ) = A*B* C DB* ( -1 3 • /70** ( BA3+ AB3 ) + * 02/ ATS ) + 1 69 0 * 

1 ATB*DM*C**2/( l44.*l ?25. ) ) 

XK0(20.3)=XK0( 14,9)=-XK0(21 »2 ) 

XKO (?1 ,3 )=XKO (15»9)= DB* C 54 • /35* *B A3- 1 56 • /35« * A83-72. / < 25 .* ATB ) > 

1 +1 1 7.*ATB*DM*C**2/(?.*1225) 

XK0<?2*3)=XK0<21 .4 > = A* ( OB* ( 27* /35 • *BA3-22 ./35 .* AB 3- 1 ./ATB* ( 6 */25* 
1 +6**XMU/5* > >+33.*ATp*DM*C**2/(4**l 225* > ) 

XKO (23 1 3 ) =XKO, ( l 7 * 9 ) = B* ( DB* ( - 1 3 • /35 • *B A3+78 . /35 * * AB3+6 */ C 25 a *ATB ) ) 
1 -169»*ATQ*DM*C**2/< 12.*! 225. > ) 

XKO (24. 3 )=XKO ( 15* 12 > = A*B* ( DB* C - 1 3* /70* *BA3+1 1 ./35**AB3+1 o / ATB* 

1 (.02+. I*XMU> >-143.*ATB*DM*e**2/<72.*1225. > > 

XKO (4 .4 )=XKO (10. 1 0 )=XKO( 16. 1 6 ) =XKO ( 22 . 22 ) a A **2* ( DB* ( 52 . / 35 * *BA3 
1 +4./35.*AB3+8./(25.*ATB > ) + 1 3 . * ATB*DM*C** 2/ < 3 • * 1 225 c > ) 

XKO (5*4 ) = XKO (17*16)= A*B* (DB* C 1 l • /35. * ( BA3+AB3 ) + l a /AT8* < o 02+ 1 .2* 

1 XMU) 1+121 .*ATB*DM*C**2/(36.*1225. ) ) 

XKO ( 6*4 ) =XKO( 12. 1 0) = A **2*8* (DB* ( 22* /l OS . *BA3+2./35 .* AB3+2 • /ATB* 

1 ( 1 ./75.+XMU/15* ) )+l 1 .*ATB*DM*C**2/( 18. *1225. ) ) 

XKO (8*4 )=XK0(22.B ) = XKO (20. 1 0 ) = XK0 ( 20 . I 6 > = -X<0 UO *2 ) 

XK 0(9*4 ) = XK 0(22* 15)=-XK0( 10.3) 

XKO ( t 0.4 )=XKO (22 . 16 )= A**2* ( DB* ( 26 . /35**BA3-3 • /35* *AB3“2 o / (25 * * ATB 
1 > )-l3.*ATB*DM*C**2/(4.*1225. ) ) 

XK 0 ( l i * 4 ) =XKO (23. 16 ) = A*B*(DB*(-1 1 . /35* *B A3+ 1 3 ./70 - *AB3- 1 ./ATB* ( 

1 .02+. 1*XMU) )+143.*ATB*DM*C**2/(72.*1225. ) ) 

XKO( 12.4 )=XKO ( 1 0.6) = A**2*B*(DB* (1 1 ./I 05 • *BA 3- 3 • /70 • *AB3“ 1 e /ATB* 

1 ( 1./150.+XMU/30, ) )-i 1 . *ATB*DM*C**2/ < 24 .*1225. ) ) 

XKO () 5.4 >=XKO (22.9) =-XKO (16.3) 

XKO ( 16. 4 )=XK0 (22. 10)= A**2* ( DB* ( 9 ./ 35 • *BA 3+3 */35 • * AB3+2 . / ( 25.*ATB ) 
1 )-9.*ATB*DM*C**2/(B.*1225. ) > 

XK 0 ( i 7. 4 ) =XK 0 ( I 6 . 5 ) = A*B*(DB*( l 3 • /70 • * ( BA3+ AB3 ) - . 02/ATB ) - 1 69 • *A T0 
; 1 *DM*C**2/( 144. *1225. ) ) 

XKO( 18.4 )=XKO (24. 1 0 )= A**2*B* (DB* ( - 1 3 • /2 l 0 . *BA3- 3 . /70 • * AB3- 1 • / 

1 ( 150.*ATB ) >+l 3.*ATB*DM*C**2/(48.* 1225. ) ) 

XK0(20.4)=XK0U0.8)=XK0( 1 6 . 8 ) =XKQ ( 22 . 20 ) =“XK 0 ( 4 . 2 ) 

XKO (22. 4 )=XKO (16.10)= A**2* ( DB* ( 1 8 « /35o *BA3-4 * /35* * A83-8 • / ( 25 ** 

1 ATB ) ) + 3.*ATB*DM*C**2/<2.*i225. ) ) 

XKO (23. 4 )=XKO (16.11 >= A*B* ( DB* ( — 13. /70 . *B A3+ 1 1 ./35.*AB3+1 ./ATB* 

1 ( .02+* !*XMU> )-l 43.*ATB*DM*C**2/ (72. *1225. ) ) 

XKO (24 »4)=XK0 ( l8, 10)= A**2*B*(DB* (-1 3. /I 05 * *BA3+2 o /35<s *AB3 + 2 c / ( 

1 75. *ATB ) )-l 3.*ATB*DM*C**2/( 36**1 225* ) ) 

XK0(5*5)=XK0< 1 1 . I l )=XK0< 17* 1 7 ) = XKO ( 2 3 . 23 > = B**2* ( DB* ( 4 • /35 » *B A3 
1 +52./35.*AB3+S./(25.*ATB> ) + l 3 * * ATB*DM*C**2/ ( 3 .* I 225* ) ) 

XKO (fi. 5 >=XKO (24.23)= A*B**2* ( DB* ( 2 « /35* *B A3+22 * / 1 05 •* AB3+2 * /ATB* 

1 ( 1 ./75. + XMU/15. ) ) + l 1 o*ATB*DM*C**2/ < 18**1225* > ) 
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XKO< 7*5 )sXKO < I 1 , 7 >=XkO(23* J 9}=XK0( 19* 17 )=-XKO < 5 , J ) 

XKO (10*5) =XKO ( 22 *17) s-X<0 ( 1 1 *4) 

XKO( i 1 ,5)=XK0(23*17)= 6**2+ < DB* ( -4 - /35 . *BA3+ 1 8 ,/35 , *AB3-8. / ( 25. * 

1 ATB > )+3.*ATB*DM*C**?/(2.*l225, > ) 

XKO( 12,5)=XK0(23, 18>= A*B**2* ( DB* (2 , /35 . *BA3- 1 3* / 1 05. * AB3+2, / ( 

1 75#* ATB) )-l 3,*ATB*DM*C**2/(36,*1225, > ) 

XKO ( 1 5, 5 )=XK0 (21,11) = -*XKO (17*3) 

XKO < l 7, 5 ) s=XKO (23, 1 l ) = B**P* ( D8* ( 3 • /35 • *BA3+9, /35, # A63 + 2 * / (25,*AT3 ) 
1 )-9.*ATB*0M*C**2/(8.*l 225. ) > 

XKO ( 1 8 * 5 ) 3 XKO (23, 12)= A*B**2* ( DB* ( -3 . /70 • +BA.3- 1 3 * /Z 1 0 . *AB3- 1 ,/ 

1 ( 150**ATB ) ) + l 3«*ATB+DM*C**2/(48,*l 225, > ) 

XKO<19,5) = XKO ( 1 7, 7 ) = XKO (23,7 > = XKO ( 19,11 )=-XK0(17, 1 ) 

XK 0(^1 , 5 ) = XK0 (15,11 ) = -XKO (23,3) 

XK0(22,5)=XK0 ( 1 7, 10)=-XKO(23,4) 


XKO (23,5) -XKO (17,11 ) = 8**2+ ( DB* (-3./35.+BA3+26./35, *AB3-2,/ (25, * 

1 ATB) ) — 1 3, *A TB*DM+C**2/ (4**1 225, ) ) 

XKO (?4 ,5 ) =XKO (23,6 1= A*B**2* ( 06* < ~3* /70 . *B A3 + i 1 ,/J 05**AB3-1 • /ATB* 

1 ( 1 ,/l 50.+XMU/30* > )-l 1 ,*ATB*DM*C+*2/ <24, *1 225. ) ) 

XKO (6,6 ) = XKO ( 12,12)= XKO (18* 1 8 ) =XKO ( 24 , 24 ) = A **2*B + * 2* ( OB* ( 4 . / l 05 • * ( 
1 BA3+AB3 >+8,/(225**AtB) > + AT8*DM*C**2/ ( 9, * 1 225 , ) ) 

XK0(9,6 )=XKO (21 ♦ !8)=-XK0( 12, 3) 

XKO ()1,6) =XKO (24 , 1 7 ) a-XKO ( 1 2 ,5 ) 

XKO( j 2,6 )=XK0 (24,18)= A**2*9**2* (DB* (2,/ 1 05,*BA3-1 ,/35«* AB3-2,/ 

1 (225,+ATB ) )-ATB*DM*c**2/< 12,*1225. ) ) 

XKO ( i 6, 6 ) =XKO (22 , 1 2 ) =-XKO (18,4) 

XKO( 17,6)=XK0(24. 1 1 )=-XKO( 18,5) 

XKO ( 1 8 , 6 > 3 XKO (24,12)= A** 2 + B **2* < DB* ( - 1 • /70 • * ( BA 3 + AB3 > + 1 . / ( 450 . * 

1 ATB ) )+ ATB*DM*C**2/ (16**1225,)) 

XKO (21 , 6 ) = XK 0 ( 1 8 , 9 > = — XKO (24,3) 


XKO <5>2,6)=XK0 ( 16, 12) = - XKO (24,4) 


XKO (24,6 >=XKO (18,12)= A**2*B**2* ( DB * ( - 1 • /35 • *B A3+2 , / 1 05 , * AB3-2 . / 
1 (22s,*ATB))-ATB*DM*C**2/( 12, *1225. ) ) 

XK 0 ( l 0 ♦ 9 ) = XK 0(16,15)= — XK 0(4,3) 

XK0< 12. 9)=XK0 (24,21 )=-XK0(6,3) 

XKO ( 16,9) =XK0 (15,10) =— XKO (22,3), 

XK 0 ( ?4 ♦ 9 ) = XK 0(21 , 1 2 ) =-XKO ( 1 8 , 3 ) 

XKO ( l 1 , !0>=XKO<23,22)=-XKO(5,4 ) 

XK0(?3, 10 ) = XKO (22, l 1 ) = -XKO( 17,4) 

XKO ( 12, 11 ) = XKO ( 18*17) = — XKO (6,5) 

XKO ( 18, 11 ) = XK0(17,12) = — XKO (24,5) 

XKO ( 15, 14 )=XK0(21 1 20 ) = -XKO ( 3 ,2 ) 

XKO( 17, 15 )=XK0(23,21 >=-XK0(5,3) 

XKO(pO, 15)=XK0(21 ,14 >=-XK0(9,2 ) 

XK0( ? 3, 15) = XKO ( 2 1 ,17)=-XK0(1 l +3) 
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XKO { 1 8« 16 ) = XK0(24 *22 )=-XKO <6 *4 ) 

XK0{?4« 16)=XK0(22*l8) --XKO ( 12*4) 

C MAKE USE OF SYMMETRY TO COMPLETE THE MATRIX 

DO 5 J= 1 *24 « 

K = J ♦ 

DO 5 I=K,24 * 

I F ( ! *EQ * J ) GO TO 5 ♦ 

XKO ( J* I > =XKO ( I * J ) « 

5 CONTINUE . 
return • 

END ♦ 

SUBROUTINE REDFORC ( JoOF *NE ) ♦ 

COMMON /ZZ5/JR (150) 

5 /ZZ6/F0RC (21 6* 1 1 

C THIS SUBROUTINE ELIMINATES FORCES AT REACTIONS 

DO 4 3 I N= 1 • NE 
LL=JR( I N ) — I N+ I 
LEF = jDOF— 1 
DO 43 N=LL*LEF 
FORC <N« 1 ) =FORC ( N+ 1 * 1 >♦ 

A3 CONTINUE 
RETURN 
END 

SUBROUTINE ADDUP (MN,XMK» XKO ) 

C0MM0N/ZZ3/N1 ( 50 > *N2 ( 50 ) *N3 ( 50 > » N4 ( SO ) 
DIMENSION XMK(216«216> *XK0<24»24 ) 

LI =N1 (MN) « 

L2=N?(MN) « 

L3=N3CMN) n 
L4 = N4 ( MN ) * 

DO 4 1=1*4 , 

DO 4 J= 1 « 4 « 

N I =6*L 1 —5 , 

NJ=6*L ! -5 . * 

I F ( I ,EQ • 2 > N!=6*L2-5 , 

I F ( J ,EQ • 2 ) NJ-6*L2— 5 ♦ 

I F ( I *EQ • 3 ) NI=6*L3-5 
I F ( J*EQ • 3 ) NJ=6*L3-5 « 

IFM.E0«4) N I =6*L4— 5 « 

IF(J.EQ.4) NJ=6*L4-5 « 

DO 4 K I =1 ,6 « 

DO 4 KJ=1 * 6 « 

KKI=nI4-KI-1 , 

KKJ=NU+KJ-1 . 
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KXI=<S*!-6+Kt * 

KXJ=6*J-6+KJ * 

XMK(<KI *KKJ ) =XMK (KK I « KK J ) + XKO ( KX I *KXJ) 

IF ( (K<I -EQ.27).AND» (kKJ.EQ.28) > 2*4 

2 PRINT 3*XMK(KKI *KKJ) 

3 FORMAT (/40X*E20.4/) 

4 CONTINUE ♦ 

RETURN * 

END « 

SUBROUTINE W'ASH ( JDOF * NE » XMK ) 

COMMON/ZZ5/JR U 50 ) 

DIMENSION XMK (21 6 *21 6) 

C SET DESIRED ROWS AND COLS TO ZERO 

DO 31 N= 1 *NE « 

LL=JR(N> « 

DO 3t 1=1 « JDOF * 

XMK (LL * I ) =0 • 0 * 

XMK ft *LL)=0.0 • 

31 CONTINUE * 

C COMPACT MATRIX ♦ 

DO 33 I N= 1 * NE « 

LL=JDf IN1-IN+! « 

LEF= JDOF— 1 « 

DO 33 N=LL*LEF * 

DO 32 1=1 ♦ JDOF . 

XMK ( N ♦ I 1 = XMK ( N+ 1*1) 

32 CONTINUE 
DO 33 1=1* JDOF 
XMK < J « N > - XMK ( I « N+ 1 ) 

33 CONTINUE 
RETURN 
END 

SUBROUTINE EXPDEF ( JDOF * NE > « 

COMMON/ZZ5/JR (150) 

7 /ZZ7/X (216) 

THE PURPOSE OF THIS SUBROUTINE IS TO ARRANGE THE DISPLACEMENTS Ih- 
THE ORDER OF THE NODE NUMBERING 
DO 201 I N = 1 * NE 
LF - JOOF— JR ( IN) 

DO 201 N=1 *lf 
X ( JDOF-N+ 1 ) = X ( JDOF—N ) 

201 CONTINUE 

DO 2r>3 11 = 1 *NE 
LL= JR ( II)* 
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X (LL ) =0 « 0 « 

CONTINUE 

RETURN 

END 

SUBROUTINE STRESS (Mm). 

COMHON/ZZI /AZ (30) , 9Z ( 30 ) * CZ ( 30 ) * DMZ ( 30 ) . DBZ < 30 ) *XMUZ(30) 

2 /ZZ3/N1 (50>.N2(50) .N3 ( 50 ) ♦ N4 ( 50 > 

6 /ZZ7/X ( 2 1 6 ) 

7 /ZZQ/S <3.24 ) . 

a /ZZ9/STRSR (3 ) * 

9 /ZZ1I/XE(24) 

9 /ZZ1 2/BTA1 .3TA? 

THIS SUBROUTINE COMPUTES THE STRESS MATR I X « S ♦ THE MULTIPLIES IT. 
BY THE DISPLACEMENT VECTOR FOR THE ELEMENT.. 

COMPUTE EACH TERM OF THE ELEMENT STRESS MATRIX. 

A = AZ (MN ) 

B = BZ (MN ) 

C=CZ(MN) 

DM=OmZ ( MN ) 

DB = DF»Z (MN) 

XMU=XMUZ (MN ) 

Sll.t ) =S ( 1 ♦ 1 9 ) *-DM/ (?.*A ) * 

S ( 1 . 2 ) =S ( 1 . 8)=-XMU*DM/(2.+B> . 

S (1 .3 ) =S C I .9 ) =S ( 1 . 1 5 ) =S ( 1 .21 > = XMU*DM*C/4. . 

S ( I .4 )=£( 1 . 22 )=A*S <1 . 3 )/4. . 

S(1.5)=S<1*11 )=B*S(1 . 3 ) /4 . f 
S(!»6)=S(lt| 8 )=A*B*S ( 1 .3 >/l 6.. 

S( 1 . 7 ) = S ( 1 . 1 3 ) =— S ( 1 . 1 ) . 

S( 1 ♦ 10 )=S< 1 ♦ 1 6 ) -— S (1.4), 

S ( 1 » 1 2 ) = S < 1 . 24 ) =-S ( 1 .6). 

S ( 1 « t 4 ) = S ( 1 . 20 ) =-S ( 1 .2). 

S ( 1 . 17) =S(1 « 23 ) =— S ( 1.5). 

S ( 2 . 1 ) = S (2. 1 9)=XMU*S (1.1 ) » 

S (2.2 )=S (2. 8) --DM/ (2**8 ) . 

S(2.3)=S(2,9)=S(2.15)=S(2«2! >*DM*C/4. « 

S ( 2 ♦ 4 >=S(2.22)=A*S(2.3)/4.. 

S(2*B)=S(2« 1 1 )=B*S(2.3>/4.. 

S (2.6 )=S(2* 18)=A*B*S (2.3)/16.. 

S( 2. 7 > = S( 2* 13) s-S (2.1 > * 

S<2. 10)=S(2.I6)=-S(2.4) . 

S (2. t 2 ) =S (2.24 )=-S{2.6 ) . 

S(2.14)=S(2.20)=-S(2.2), 

S (2. I 7) *S (2*23 ) «-S(2 ,5 ) . 

S ( 3 . | ) =S ( 3 . 7 ) =— ( 1 .-XMU)*DM/ (4**B >. 
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S<3,2)=S<3*20 )=-( I • — XMU )*DM/(4**A ) , 

S (3*8 ) =5 ( 3* 1 4)=-S(3.2 ) * 

S(3« 13>=S(3*19)=~S(3«1 >. 

S(3,3)=S(3,4>=S(3,5>=S(3,6)=S(3*9)=S<3.10)sS(3,ll>=S<3,12)=S{3,15), 
1 =S(3,1 6)=S(3» 17)=S<3, 18 >=S(3,2l )*S(3,22 ) =S < 3 ♦ 23 ) =S < 3 , 24 )=0*0* 

C PRINT out element stress matrix* 

PRINT' 751 * MN, 

751 FORMAT C40X* ^ELEMENT STRESS MATRIX FOR ELEMENT NO *,13/), 

PRINT 753* ( (S(I»J)*J=1 i24) , 1 = 1*3)* 

753 FORMAT ( 6E20 .4 ) ♦ 

CONSTRUCT THE DISPLACEMENT VECTOR FOR THIS ELEMENT* 

DO 700 ! E= 1 *24 , 

I F U* 7 . GE * 7 ) GO TO 70?* 

ND!=ft*(NI (MN)-l >« 

XE< I F ) = X ( ND 1 -f I E ) * 

GO To 700* 

703 IFUf.GE. 13) GO TO 7o5 * 

ND2-6* < N2 ( MN ) — 1 ) , 

XE ( IF }=X(ND2+IE-6 ) « 

GO TO 700. 

705 IFUf.GE. 19) GO TO 7o7* 

N03=6*CN3(MN>-I ) 

XE ( IF ) =X ( ND3 + I E— 1 2 ) « 

GO TO 700. 

707 ND4 = fi* ( N4 < MN ) — 1 )« 

XE< I F ) = X < ND 4+IE— 18) . 

700 CONTINUE, 

PRINT 709, MN, 

709 FORMAT </40X,*D I SPLACFMENT VECTOR FOR ELEMENT NO *,I3/), 

PRINT 71 1 « ( XE ( T > , 1=1 ,24 > « 

711 FORMAT (4X.F20.4), 

BTM = (XE(3)-XE<9)+XE<21 >-XE(15> )/(2*A> 

BTA2=(XE(3 ) — XE ( 2 1 )+XF<9)-XE(!5) )/(2.*B)+C*(XE(2)+XE<8 )+XE ( 1 4 )+XE ( 2 
1 0 ) )/4. 

PRINT 720 

720 FORMAT < /40X , *A VERAGE ROTATIONS ABOUT Y AND X AXES*/) 

PRINT 722.3TA1 ,BTA2 
722 FORMAT ( 40X , E20 .4 ) 

MULTIPLY STRESS MATRIX BY THE DISPLACEMENT VECTOR TO OBTAN THE. 
INPLANE STRESS RESULTANTS , ENX * ENY AND ENXY, 

CALL MATRIX(20,3,24*1 , S , 3 , XE , 24 ♦ STRSR, 3 ) , 

PRINT 713.MN* 

713 FORMAT f /40X**STRESS RESULTANTS FOR ELEMENT NO*. 13/), 

PRINT 715. (I, STRSR (T)*I»1*3>* 



715 FORMAT (4X.*STRSR<** T 3 * * ) =* « El 2 • 4 ) . 

RETURN* 

FND * 

SUBROUTINE ELEMK1 ( MN > ♦ 

COMMON/ZZ l/AZ (30).BZ(30)*CZ<30> «DMZ ( 30 > *DBZ ( 30 ) » XMUZ t 30 ) 

8 /ZZ9/STRSR ( 3 ) * 

9 /ZZ1 0/XK1 (24*20 > « 

A=AZ(MN) 

B=BZ(MN> 

C=CZ(MN> 

OM=DMZ ( MN ) 

OR=OPZ ( MN ) 

XMU=yMUZ(MN> 

C ZERO OUT XK I « 

DO 900 1=1 .24* 

DO 900 J= 1 *24 * 

XX 1 ( I * J 1=0.0* 

900 CONTINUE* 

SI =STRSR( 1 1 * 

S? = ATRSR (2 ) * 

SI 2 = STRSR ( 3 1 * 

XK 1 (2*2) = XK 1 (8*0) =XKl (14*14 )=XK1 ( 20 * 20 ) =S2*A*B*C+*2/9. * 

XK 1 (3*3) =XK 1 ( 1 5* 1 5 ) *S1 *8/ (3«*A)+S12/2*+S2*A/(3**B) ♦ 

XK 1(9*9)= X< 1 (21*21 ) = <:l*B/<3**A)-S12/2* + S2*A/(3**8)* 

XK 1 (3*2)= XK! (2*3)=(Sl 2*B-f S2* A ) *C/6 • • 

XK! (8*2) =XK 1 (2.8 > =S2 *A*B*C**2/ 1 B . * 

XK 1 (9*2) =XK 1 ( 2 . 9 ) =- ( 5 1 2*9-S2*A /2 • >*C/6* . 

XK1 ( 1 4*2) =XK 1 (2 * 14 ) = X<1 (2. 2 )/4. « 

XK 1 ( )5*2)=XK1 (2*1*5) = -XK1 <3*2)/2. ♦ 

XK 1 ( pO ♦ 2 ) =XK 1 (2*20)=XK1 ( 14 *8)=XK1 (0* 14 )=XK1 (20*14) = XK 1 (14*20) = * 
1 XK) (3.2 )/2. * 

XK 1 ( 21 * 2 )=XK 1 (2*21 ) = ( S 1 2*B/2 .-S2*A )*C/6* » 

XK 1 (R.3)=XK! (3*B)=(Sl2*R+S2*A/2. >*C/6.* 

XK1 (9*3) =XK 1 ( 3* 9 > =-Sl *B/ < 3* * A ) +S2*A/ (6**B) * 

XK 1 ( ! 4 * 3 ) =XK1 ( 3 ♦ 1 4 ) = -X«l ( l B * 2 ) * 

XK 1 ( ! 5* 3 ) =XK1 (3* 15 ) =-Sl*B/ ( 6**A )-Sl2/2«-S2*A/(6**B ) • 

XK 1 ( pO ♦ 3 ) =XK 1 C3»20)=(S12*B/2*+S2*A )*C/6*. 

XK1 (31 * 3 ) =XK 1 (3.21 ) = 51 *B/ ( 6 • * A ) - S2*A/ ( 3 • *8 > * 

XK! (9*9) =XK 1 (0.9 )=- <S! 2*B-S2*A )*C/6. ♦ 

XK 1 ( 1 5 * 8 ) =XK 1 (8* 15>=-(S12/2.*B+S2*A)*C/6. * 

XK 1 ( 30 * 8 )= XK 1 <8»20) = xKl(2.2)/4.* 

XK1 (?1 *0)=XK1 (8*21 ) = ( S 1 2*B— S 2*A ) *C/ 1 2 • * 

XK l ( 1 4 » 9 ) =XK 1 (9*14) =— XK 1 (21*2)« 

XK 1 ( l5.9)=XKl (9*15)=X<1 (21 *3>* 
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XKl (?0, 9>=XK! <9«20>=XK1 <9*8)/2* * 

XK1 (21 *9)sXKl (9*21 ) =-S 1 *B/ < 6 * * A ) +S l 2/2 • -S2*A/ ( 6 • *B ) « 

XKl (20* 15>=XK1 (1.5*20 ) =- ( S J 2*B+S2*A/2 • >*C/6« * 

XKl (?! ♦ 14)= XKl (14*21 )=-XKl (9.21* 

XK 1 < t 5* 14 )=XK 1 (14*15)= — XKl (3*2>* 

XKl (21* 15) =XK 1 (15*21 )=-Sl*B/(3**A)+S2*A/(6**B>* 

XKl (?1 * 20 ) = XK 1 (20*21 > = -XKl (3*2)* 

RETURN* 

END* 

SUBROUTINE' ELEMK2 ( «N ) 

C0MM0N/ZZ1 /AZ ( 30 ) *SZ (30 ) *CZ ( 30 ) ♦ DMZ < 30 ) «D6Z ( 30 ) *XMUZ < 30 ) 

9 /ZZ-1 2/BTA1 .BTAp 

Q / ZZ 1 3/XK2 (24*24) 

9 /ZZ 1 4 /SO 

A=AZ(MN) 

B=BZ ( MN ) 

C=CZ(MN) 

DM=OMZ(MN) 

OB=DRZ(MN) 

XMU=yMUZ ( MN ) 

C ZERO OUT XK2 ( I * J) 

DO 9*0 1=1.24 

DO 950 J=1 *24 
XK2(I*J)=0*0 
950 CONTINUE 
C DEFINE constants 

Cl =1 4 o-xmu 

C2=1*0+XMU 
C3= 1 # 0 — 3 * 0*XMU 

c define second order terms* a modulating factor, so* will determine 

c wheather the second ORDER terms will be used or not 

ALF 1 sSO* OTA 1**2 + C1 *RTA2**2/2. 0 ) 

ALF2=S0* (BTA2**2+C1 *pTAl **2/2*0 ) 

ALF3=S0*C2*BTA1*BTA2 

XK 2(9*1 )=XK2( 1 *2 >=-C#< •S*C1*A*BTA1+xMU*B*BTA 2)/6.0 
XK2 ( R * 1 )=XK2 ( l ,3)=-( .5*C1*A1 B+B/A)*BTA1 /3*-C2*BTA2/8*0 
XK2 (8*1 )=XK2(1 *8 )=-C* ( .25*01 *A*BTA1+XMU*8*BT A2 )/6*0 
XK 2(9*1 )=XK2(1 ♦ 9 ) =- ( , 25*C 1 * A /B-B/ A ) *BT A 1 /3 . 0+C3*BT A2/8 • O 
XK 2(14*1 ) = XK 2 < 1 * 14 )=*5*XK2 (2*1 ) 

XK2 ( 1 5* 1 ) =XK 2 (1 « 1 5 ) = ( *5*C1 *A /B+B/ A ) +BTA 1 /6*0+C2*BT A2/8* 

XK2 ( e>0 * 1 >=XK2(1 ♦ 20 ) = —C* (Cl *A*BT A 1 +XMU*B*BTA2 ) / 1 2 • 

X<2(21 * 1 ) = XK 2(1 *21 ) = ( C 1 *A/B-B/A ) *BTA 1 /6. -C3*BTA2/8 • • 

XK2(?*2 )=-C*( *5*C1 *B*BTAl+A*BTA2 )/3.+A*B*C**2*ALF2/9* 

XK2 <3*2 ) =XK2 (2*3 )=-C2*BTAl/S*-( • 5*C 1 *B/A+ A/0 > *BTA2/3 • + C* < *5*B*ALF3 
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1 +A*ALF2 )/6#+C**2*A*B*BTA2/9. 

XK 2(7*2 ) =XK2 (2.7) =-C* ( • 25*C 1 *A*BTA 1 -XMU*B*BT A2 )/6* 

XK2(8,2 >=X<2 (2.9)=-A*0TAl *C/6 • +A*B*C**2* Al_F2/ 1 8* 

XK2 ( 9 » 2 >=XK2 (2.9 ) =-C.?*BTA 1 /8 • + (Cl *B/A-A/B > *BT A2/6 • -C* < B* ALF3-A* 

1 ALF2 >/! 2.+A*B*BTA2*C**2/18. 

XK 2(13,2)- XK 2(2. I 3 > =c* ( • 5*C 1 *A*BT A l -f XMU*B*BT AP ) / 1 2 • 

XK?( 1 4,2 > = XK2(2 , 14 )=a*B*ALF2*C**2/36. 

XK2( 1 5. ’2 )=XK2 (2 ♦ 1 5 ) = c2*BTA I /8 • + ( • 5*C 1 *8/A+ A/B )*BTA2/6. -C* ( .5*6* 

1 ALF3+A*ALF2 )/l 2.+A*B*BTA2*C**2/36. 

XK2( 1 9,2 )=XK2 (2 ♦ 1 9)~r* (Cl *A*BTAl -XMU*B*BTA2 ) /t 2* 

XK2(?0,2 )=XK2 (2.20)=-C 1 *B#BT A ! *C/ 1 2 • +5* *A*B* ALF2*C**2/36« 

XK2C21 . 2 )=XK2 (2.21 ) =C2*BTA 1 /8. - ( C 1 *3/ < 4 * *A ) - A/B ) *BTA2/3. + C* < « 25*B* 
1 ALF3-A4ALF2 > /6 • + 5 • * A *B*BTA2*C**2/ 36 • 

XK2 ( ? * 3 ) = C* (XMU*B*BTAl+A*BTA2 ) /3 . +ALF3/4 • +B*ALF 1 /( A»3* )+A*ALF2/ (3* 
13.) 

XK2(7.3)=XK2 (3,7 >=-(ct *A/(4,*B)-B/A ) *BT A 1 /3 . ~C3*BT A2/Q • 

XK2 ( fl. 3 > =XK2 ( 3 , 8 )s=C3*BTAl /8*+ ( C 1 *B/ A- A/B ) *BT A 2/6 . + C* ( B*ALF3+A*ALF2 
1 )/l 2 , + A*R*BTA2*C**2/l 8. 

XK2(9*3)=XK2(3,9 5 =A*rTA 2*C/6 *-B*ALFl / C 3. *A >+A*ALF2/ ( 6* *B ) 

XK2 ( i 3, 3 >=XK2 (3, 13 ) =XK2( I 5* 1 > 

XK 2 ( 1 4, 3 )=XK2 (3,14)=C2*BTAl/8*+( . 5*C 1 *8/ A+A/B >*BTA2/6*+C* < .5*B* 

1 ALF3 )/l 2.+A*B*BTA2*C**2/36* 

XK2 ( 1 5, 3 )=XK2 (3. I B ) = -ALF3/4 . “B*ALF 1 / ( 6 •* A ) -A* ALF2/ ( 6 • *B ) 

XK2 M 9 , 3 ) -XK2 ( 3 . 1 9 ) = < C 1 *A/B-B/A ) *STA 1 /6 • + C3*BTA2/8 . 

XK2(20,3 )=XK2<3, 1 9 > = -C3*BT A 1 /8 « - < . 25*C 1 *8/A >*BTA2/3.+C* ( * 25*6* 

1 ALF3+A*ALF2 ) /6 . +5 • *A*B*BTA2*C**2/36 • 

XK 2 ( 2 1 « 3 ) =XK2 (3,21 ) = xMU*8*BT A 1 *C/6 . +B*ALF 1 / ( 6 * *A ) - A*ALF2/ ( 3 • *B > 

XK2 (8,7 )=XK2(7,8 ) --C* ( Cl *A*BTA ! /2 . -XMU*B*BtA2 )/6. 

XK2 (9,7) =XK2 (7,9 ) =- ((*2* A/ (2**B )+B/A ) *BT A 1 /3 • +C2*BTA2/8 • 
XK2(14,7)=XK2(7. 14)=-XK2( 19,2) 

XK2 < 1 5 , 7 ) =XK2 (7.15) = XK2 (19.3) 

XK2(?0.7)=XK2 (7,20) =XK2( 8,7) /2. 

XK2 (21 «7 )=XK2 (7.21 ) =XK2 (9,7) 

XK2(8,9 >=C* CCl*B*BTAl/2.-A*BTA2 ) /3«+A*B*ALF2*C**2/9. 

XK2 (o, 8 )=XK2 (8,9 )=C2*BT A 1/8.- ( .5*Cl*B/A+ A/B > *BTA2/3*-C* ( *5*B*ALF3 
1 -A*AlF2 )/6.+A*B*BTA2*C**2/9. 

XK 2 ( 1 3 , 8 ) =XK 2(8*13) =->XK2 (20,1) 

XK2 (14,8) =XK2 (6. 14 ) =C 1 *B*BTA 1 *C/ 1 2 • + A*B* ALF2*C**2/ 1 8* 

XK2 ( 15.8 ) =XK2(8, 15) =-C3*BTAl /8.. - < ♦ 25*C 1 *B/A- A/B )*BTA2/3.-C* (B*ALF3 
1 /4 ,+A*ALF2 >/6 ,+A*B*RtA2*C**2/18* 

XK2( 1 9,8) =XK2( 8.19)=-. 5*XK2< 8, 7) 

XK2 (20. 8 )sXK? (8, 20)=x<2 (8,2 ) 

XK2( 21 .8 >=XK2 (8,21 ) =~C2*BTA 1 /8 • +< • 5*C 1*B/A+ A/B >*BTA2/6*+C* ( • 5*6 
l*ALF3/2«-A*ALF2 >/! 2 .+A*B*BTA2*C**2/36. 
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XK2<9»9 )=-C* (XMU*B*BTA1-A*8TA2)/3*+ALF3/4*+B*ALF1/( 3*#A )+A*ALF2/( 

1 3 • *8 ) 

XK2< i 3.9) =XK2(9« 13 1 =X*2< 2 1 • l ) 

XK2 C 1 4 « 9 ) =XK2 (9*14) =C3*BTA 1 /8.- ( *25*C1 *8/A-A /B >*BTA2/3.-C* ( • 25*B* 

I ALF3— A* ALF2 ) /6 •+A*B*pT A2*C**2/1 8* 

XK2( 15.9) =XK2(9« 1 5 > =-XMU*B*BTA 1 *C/6 •+B#ALF 1 / ( 6 .*A > -A*ALF2/ ( 3 • *B ) 
XK2( 19* 9 )*XK2 (9* 19) =X<2(2l *7) 

XK2(?0*9 ) =XK2(9*20)*-C2*BTAl/8* + ( .5*Cl*B/A+A/8 >.*BTA2/G*-C*<0*ALF3/ 
1 2.~A*ALF2 )/l2.+A*B*BTA2*C**2/36* 

XK2<?1 * 9 ) = XK 2(9*21 ) « ALF3/4 .-B* ALF 1 / ( 6* *A >-A* ALF2/ { 6 .*B ) 

XK2( 14* 13 )=XK2( I 3 * 1 4 ) = -.5*XK2 < 1 4., 1 ) 

XK2 ( t 5 * 13) -XK2 (13*15 )*XK2 C 3* 1 ) 

XK2 ( ?0i 1 3 ) =XK2 ( 13*20 )=-XK2 (8*1 ) 

XK2(?1 « l 3 ) -XK2 ( 13*21 )=XK2(9* 1 > 

XK2( 14. 14 )=C*< *5*C1 *B*BTA 1 +A*BT A 2 )/3 •+A*B*ALF2*C**2/9# 

XK2( 15*14 )=XK2( 14*15 )=-C2*BTAl/8*-( *5*0 *B/A+A/B)*BTA2/3.-C*( ,5*B* 
1 ALF3+A* ALF2 >/6 •+A*B*RTA2*C**2/9« 

XK2(f9* 14)=XK2( 14*19)=“XK2 (7.2) 

XK2(pO* 14 ) = XK2( 14.20 > = A*BTA2*C/6.+A*B*ALF2*C**2/18. 

XK2<?1 * 14 )=XK2( 14.21 ) = -C3*BTA 1 /8 . + {C 1 *B/ A- A/B ) *BTA2/6 . +C* ( B*ALF3- 
1 A#AL^2)/I2.+A*B*BTA2*C**2/18, 

XK2( 15. 15 > = -C*(XMU*B*BTAI+A*BTA2 ) /3 ♦ +ALF3/4 .+B*ALF1 / ( 3 . *A >+A*ALF2/ 
1 ( 3**P ) 

XK2(J9. 15)=XK2(15.19)=XK2(7*3) 

XK2 (20.15) =XK2 (15*20) =C3*BTA I /8 • + (C ! #B/A-A/B >*BTA2/6. -C* ( B*ALF3+A* 
1 ALF2 )/l 2.+ A#B*BTA2+C**2/1 8* 

XK2(21 ♦ 15 )=XK2( 15.21 ) =-A*BTA2*C/6.-B*ALF 1 / ( 3. *A )+A*ALF2/ ( 6. *B ) 
XK2(?0. 19)-XK2(19*20) =— XK2 (0*7) 

XK2<?1 . 19)=XK2( 19*21 )=XK2(9*7) 

XK2(20,20)=-C*C.5*Cl*B*BTAl-A*BTA2)/3.+A*B*ALF2»C**2/9. 

XK2( 21 *20 )=XK2 (20*21 )=XK2(9*8) 

XK2 (21 . 21 )sC* ( XMU*8*RT A 1 -A*BTA2 >/3* -ALF3/4 .+8* ALE 1 /< 3.*A >+A*ALF2/ 

1 (3.*R) 

DO 9fi0 1=1 * 24 
DO 9*0 J = 1 * 24 
XK2 ( I ♦ J ) =DM*XK2 ( t * J ) 

960 CONTINUE 
RETURN 
END 

OVERLAY (L1NK.2.0) . 

COMMON/221 /A 2 (30 ) *BZ(30) ♦ C2 ( 30 ) ♦ DMZ ( 30 j *DB2 (30 ) *XMU2( 30 ) 

1 /ZZ2/XK0 ( 24 *24 ) * ' 

2 /ZZ3/N1 (50 ) *N2 (5.0) *N3( 50 ) *N4 ( 50 ) 

4 /2Z5/ JR ( 1 50 ) - 
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5 /ZZ6/F0RC (216*1 ) 

6 /ZZ7/X (216) 

7 /ZZ8/S<3*24 )* 

8 /ZZ9/STRSR ( 3 ) » 

9 /ZZ10/XK! (24*24 > « 

9 /ZZlt/XE(24> 

9 /ZZ12/BTA1 *BTA? 

9 /ZZ13/XK2(24»24> 

9 /ZZ14/S0 

DIMENSION XMK <1 80 * 1 80 ) *XMK2( 180*180)* 

1 RTR< 1 80 > *RT ! ( 1 80 ) * lR(jN( 1 80 ) *P( 180 > *NDEX (180) 

REWIND 2* 

Read TAPE 2* ( <XMK2< I , J ) * I = 1 * 180 > • J»1 * 180 > *NO* JDOF*NE* JDOFR* 

C ADD XMK1 AND XMK2 ANn STORE IN XMK 

REWIND 3* 

READ TAPE 3 ♦ < (XMK <I*J)*I=1* JDOFR ) * J= 1 * JDOFR ) 

C MULT THE MATRIX < XMKl +XMK2 > BY XMK- INVERSE 

CALL MATRIX (20* JDOFR* JDOFR* JDOFR* XMK2 * 1 80 • XMK * 1 80* XMK • 1 80 ) 

C. GET rlGEN VALUES FOR BUKL 

NPLUF- JDOFR+1 

CALL REIG (XMK* JDOFR* 3* 3«RTR*RTI • XMK • 1 80 *NDEX * I RUN*P , NPLUS * XMK2 ) 
PRINT 171 

171 FORMAT C/45X * *E IGEN VALUES*/) 

PRINT 173* (RTRC I ) * |=l *3) 

173 FORMAT (3E20*4 ) 

RETURN* 

ENp 



